UTEC ME 87-034 


Final Report 
to the 

National Aeronautics and Space Administration 
Grant NAG-1 -283 


ACOUSTIC PROPAGATION IN A 
THERMALLY STRATIFIED ATMOSPHERE 


W. K. Van Moorhem 

Mechanical and Industrial Engineering Department , 
University of Utah 
Salt Lake City, Utah 841 1 2 


September, 1987 


l 1812351 ACCtJSTIC EECPAG1TIOB IH A 

wiSi ainnii 

total, am..) ii4 p »,aiu ms bc 
101 


G3/71 


N87-27490 

U&clas 

00928S7 



ABSTRACT 


Acoustic propagation in an atmosphere with a specific form of a temperature 
profile has been investigated by analytical means. The temperature profile used is 
representative of an actual atmospheric profile and contains three free parameters. 

Both lapse and inversion cases have been considered. Although ray solution have 
been considered the primary emphasis has been on solutions of the acoustic wave 
equation with point source where the sound speed varies with height above the ground 
corresponding to the assumed temperature profile. The method used to obtain the 
solution of the wave equation is based on Hankel transformation of the wave equation, 
approximate solution of the transformed equation for wavelength small compared to the 
scale of the temperature (or sound speed) profile, and approximate or numerical 
inversion of the Hankel transformed solution. The solution display the characteristics 
found in experimental data but extensive comparison between the models and 
experimental data has not been carried out. 
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1.0 INTRODUCTION 


Although the propagation of acoustic signals through the atmosphere has been 
studied for many years, and most atmospheric effects are understood in the qualitative 
sense, quantitative modeling of most of these effects has become an area of interest 
only recently. The dominant effect occurring in atmospheric propagation is the 
spreading of acoustic energy associated with a wave propagating in three dimensions 
over an ever increasing area, the well known spherical spreading effect, which occurs 
in an isothermal, unbounded atmosphere. In addition to this, acoustic waves are 
absorbed by the atmosphere, reflected and absorbed at the ground surface, scattered 
by turbulence, and refracted by both wind and temperature gradients. 

This report summarizes a project to develop models for the propagation of 
acoustic signals from a point source above a finite impedance ground surface in the 
presence of temperature gradients in the atmosphere. The situation of interest is the 
case of sound from a source located within a few meters of the ground propagating to a 
receiver located within a few meters of the ground through the temperature gradient 
that commonly occurs just above the ground surface. Best[1], Gieger[2], and 
Reynolds[3] all discuss the temperature gradient in this region. Within one to two 
meters of the ground the temperature generally goes through a diurnal cycle with a 
lapse condition, temperature decreasing with height, occurring in the afternoon and an 
inversion condition, temperature increasing with height, at night, see Figure 1.1. 

Shortly after sunrise and sunset the atmosphere goes through a nearly isothermal 
period when the transition from lapse to inversion or inversion to lapse condition is 
under way. This simple picture of the very complex atmospheric dynamics near the 
ground can be upset by significant winds which increase the mixing near the ground 
surface and tend to lead to a more isothermal situation, or to an overcast which can 


prevent a strong lapse condition from developing by blocking the insolation or prevent 
an inversion from occurring by blocking the radiation from the ground to the night sky. 

References 1 ,2 and 3 all discuss the classic logarithmic temperature profile given 
by 



(1.1) 

which is based on empirical results. This result, although fitting the experimental data 
well, certainly is not reasonable either for heights very near the ground or very far 
above the ground. In addition, logarithmic functions are generally more difficult to deal 
with in an analysis then are algebraic functions. For this latter reason the profile used 
in this study is 


T=T + — AT ■ 
08 1 +az 


( 1 . 2 ) 

This form of the temperature profile is shown in Figure 1 .2 along with some 
temperature data obtained by Butterworth[4]. The agreement between the data and 
the assumed function fitted to this data is excellent. Also as compared to (1 .1 ) the 

physical meaning of the parameters in (1 .2), T^, AT, and a, are clear. The assumed 
temperature profile asymptotically approaches the temperature T*. high above the 


ground. At the ground the temperature is T w + AT, thus the change in temperature 
between the ground and far about the ground is AT. The derivative of temperature 


2 



with respect to height evaluated at the ground surface is - a AT. Thus 1/a is the scale 
over which the temperature change AT occurs. For example at a height z = 1/a 


one-half of the total temperature change AT has occurred. The temperature profile of 

(1 .2) can be used to represent either a lapse or inversion condition. For a lapse 


condition the parameter AT is positive and for an inversion it is negative. 

The equation governing the wave motion is the simple acoustic wave equation 
with a sound speed varying with height and with a point source term, 

_1_ i| . V 2 p . ± e M 5( r , 5(2 . s) 

a(z) at * r 


At the ground surface, z = 0, a normal impedance boundary condition 


(1.3) 


p = -Z-|2- 
ipco 


(1.4) 


is assumed. High above the ground, z -> », only outgoing waves are permitted, a 
radiation condition. At the source height, z = s, the pressure field is to be continuous, 
and to satisfy the conditions implied by (1 .3). 

Section 2 contains a discussion of the acoustic rays that characterize both the 
lapse and inversion cases for the assumed temperature gradient, (1.2). This both 
yields a quantitative understanding of the propagation phenomena, and plays an 
integral part in understanding the modeling that follows. The model for the lapse 
condition is developed in Section 3 and the inversion model is described in Section 4. 


Section 5 contains a discussion of the conclusions developed during the project. 


2.0 ACOUSTIC RAYS 


Acoustic ray tracing is a relatively simple procedure for an axisymmetric case 
which yields a great deal of qualitative information about a given propagation situation. 
Only a brief discussion is given here. More detail is given in [5]. For an horizontally 
stratified atmosphere the acoustic rays may be determined from an acoustic form of 
Snell's law 


Cos 9(z) Cos 9(s) 
a(z) " a(s) 


( 2 . 1 ) 

where the source is located at the height s. Here 0(z) is the angle between a ray and 
the horizontal at the height z, see Figure 2.1 . Thus the right hand side of (2.1) is a 

constant for a ray emitted from the source at an initial angle 0(s). Using (1.2) to obtain 


a 2 (z) = a 2 ( 1 + 

oo 


AT 1 
T «(1+cx z) 


) 


( 2 . 2 ) 

which describes the sound speed as a function of height, the slope of a ray initially 
emitted from the source at an angle 0(s) is determined from (2.1) to be given by 


dz / Az + B' 
d7 = ± V Cz + D 

(2.3) 
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where 


A = a [ 1 + as + — - ( 1 + as) Cos 2 0(s ) ] 


B = [ 1 + CXS -K - ( 1 + Cts)(1 + 'Y”) Cos 2 e(s) ] 


(2.4) 


(2.5) 


C = a ( 1 + as ) Cos 0(s) 


( 2 . 6 ) 


and 


D = ( 1 + y -)(1 + as) Cos 2 0(s) 

oo 

(2.7) 

Equation (2.3) can be integrated to obtain the ray paths. Different results are obtained 
in the lapse and inversion cases and these will be considered separately in the next 
two sections. 

2.1 Lapse Case 


In the lapse case the quantity A is positive and integrating (2.3) to obtain the rays 



yields four different cases. For rays going upward from the source (using the positive 
sign in (2.3)) 


r-F(z, 6(s)) - F(s, 6(8)) 

for rays going downward initially from the source (using the negative sign) 

r-F(s, 0(s))-F(z, 6(8)) 


( 2 . 6 ) 


(2.7) 


for rays that were initially going downward but have been reflected upward at the 
ground (using the positive sign and (2.7)) 

r = F(z, 0(s)) + F(s, 6(8)) - 2 F(0, 6(8)) 

( 2 . 8 ) 

and for rays that initially were going downward and were refracted upward before 
reaching the ground (again using the positive sign and (2.7)) 

r = F(z, 0(s)) + F(s, 0(s)) 


The function F(z, 0(s)) is given by 


F(z, 0(s)) 

— A + 


2A/AC 


In 




(2.9) 


where 


E= AD- BC 

= a ( 1 + a s + y- ) ( ) ( 1 + a s ) Cos 2 0(s) 

oo oo 

( 2 . 11 ) 

and 


4 > = 



C( Az+B) ’ 
| A| (Cz + D) 


( 2 . 12 ) 

The absolute value of A in (2.12) is immaterial in the lapse case where A is always 
positive but significant in the case of an inversion where A may change sign. 


The rays are identified by the parameter 0(s), the initial angle at which the ray 
leaves the source. Thus a ray initially propagating downward and identified by a 


particular value of 0(s) will either be reflected upward at the ground or refracted upward 
at a turning point. In either case the reflected or reflected ray will be identified by the 


same value of 0(s) as the initial ray. Figure 2.2 is an example of the rays calculated 
from (2.6) to (2.9). 

Setting 0(z) = 0 in (2.3) yields an expression for the height at which an initially 
downward propagating ray becomes horizontal, the turning point, as 


(2.13) 




Solving this expression for Cos 0(s) yields 


Cos 0(s) = 

(2.14) 

which identifies the ray having a turning point at a height z lp . The ray that grazes the 
ground and is the boundary between reflected and refracted rays can be found by 
setting Zjp = 0 in (2.14). The ray that divides the initially upward and downward 

propagating rays is identified by 0(s) = 0. 

For the ray that grazes the ground, and is identified by the value of 0(s) defined by 
(2.1 4) with Zfp = 0, the function F(0, 0(s)) = 0 and thus this ray is defined by either (2.8) 
or (2.9). Similarly for 0(s) = 0 the ray can either be obtained from (2.6) or (2.9) since on 

this ray F(s, 0) = 0. At a turning point F(Ztp, 0(s)) = 0 when 0(s) is given by (2.14). 

The shadow boundary is more difficult to locate. It is bounded by refracted rays 
(2.9) and at a fixed height z is the maximum possible value of r for rays with turning 
points below that height. Thus the ray tangent to the shadow boundary or caustic at the 
height z is identified by solving the equation 
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a r _ a F(z, e(s)) t a F(s, e(s)) 
ae(s) ae(s) ae(s) 


(2.15) 

for 0(s) and then using that value in (2.9) to determine the location. The derivative in 

(2.15) can be obtained from 


where 


a_F 

ay 


{(1 +az) (1 -y 2 ) — 3 y 2 } 


r ip 

\/ t: 


a ( 1 - 1 2 ) 2 / 1 + a z - 'Z ( 1 + a z 4~- ) 


■^“[1 + 2 V 2 ] 

OO 

1 

2[ i -r 2 ) 2 


In 


<-Hf) 


y= 


1 + az 


T 


, AT 

1 + az+ — 


Cos 8(s) 


(2.16) 


(2.17) 

Due to the complexity of this expression an analytic solution is not possible and either 
a numerical solution of (2.15) must be obtained or the approximate relations given in 
[5] must be used. 


2.2 Inversion Case 


As in the lapse case several different types of rays occur in the inversion case. 

Integrating (2.3) with AT/T^ negative also yields three different forms for the function 

which determines the rays depending upon whether the quantity A given in (2.4) is 
positive, negative or zero. The meaning of these three cases is discussed below. For 
rays that are initially angled upward (using the positive sign in (2.3)) 


r = Fj(z, 0(s)) - Fj(s, 6(s)) 


(2.18) 

where i = 1 , 2, or 3 depending on the initial angle of the ray leaving the source. Rays 
with i =1 leave the source a sufficiently large angle upward so they do not have a 
turning point and are never refracted downward. The case i =2 corresponds to the 
limiting ray that has a turning point at infinite height. Rays described with i = 3 have 
turning points at finite height and are alternately refracted downward and reflected 
upward at the ground. These are the rays trapped by the inversion. 

Rays that initially are angled downward are given by (using the negative sign in 

(2.3)) 


r * Fj(s, 0(s))-F.(z, 6(8)) 


(2.19) 

for all three cases before they are reflected upward at the ground. The reflected 
waves are given by (using the positive sign) 


r = F.(z, 6(s)) + F.(s, 0(s)) - 2 F.(0, 0(s)) 


( 2 . 20 ) 

in all three cases. Note that after reflection the i = 3 rays are refracted downward and 
reflected upward from the ground repeatedly. In the case of these i = 3 type rays four 
more forms exist. For rays that initially were angled upward (using the negative sign) 

r = -F 3 (z, 0(8)) - F 3 (s, 0(s)) - 2 n F 3 (0, 0(s)). 

(2.21) 

after they have been refracted downward and have been reflected n times from the 
ground. For rays that initially were angled upward and have been reflected upward n 
times from the ground and have not been refracted downward following that reflection 

r = F 3 (z, 0(s)) - F 3 (s, 0(s)) -2 n F 3 (0, 0(s)) 

( 2 . 22 ) 

Thus an i = 3 type ray leaving the source upward is first described by (2.1 8) or (2.22) 
with n = 0 before it is refracted downward through a turning point. After it is refracted 
downward the first time it is described by (2.21) with n = 0. Following its first reflection 
from the ground it is given by (2.22) with n = 1 , then by (2.21 ) with n =1 between 
refraction through a turning point and reflection, then (2.22) with n =2, etc. 

For rays that initially were angled downward (using the positive sign) 


r = F 3 (z, 0(s)) + F 3 (s, 0(s)) -2 n F 3 (0, 0(s)) 


(2.23) 

after they have been reflected upward n times from the ground. For rays that were 
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f 

i 




initially angled downward (using the negative sign) 

r = - F 3 (z, 8(s)) + F 3 (s, 6(s)) -2 n F 3 (0, 0(s)) 

(2.24) 

after they have been reflected n times from the ground and have been refracted 
downward through a turning point. 

Thus an i = 3 ray that is initially angled downward at the source will first be 
described by (2.19) or (2.24) with n = 0 until it reflects from the ground, then by (2.20) or 
(2.23) with n = 1 between reflection and refraction through a turning point. Following 
the turning point and before the second reflection (2.24) with n = 1 . Then by (2.23) with 
p. = 2, etc. 

The functions F ( are given by 


F,(z. 0(8)) 


V(Az + B)(Cz + D) 

A + 


2AV AC ' $-1 ’ 


(2.25) 


for i = 1 , 


Fo(z, m = 


3CfB 


(Cz+D) 


for i = 2, and 


(2.26) 
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F a (z. 0(s)) 


J(Az + B)(Cz + D) I 

A + A / 


AC 


Tan 


-i 



(2.27) 

for i = 3. As discussed above, the i = 1 case occurs for A greater then zero or for 
Cos 0(s) < Cos © (0 > ©) where 



(2.28) 


^ Rays with values of the initial angle, 0(s), greater the value of the angle given by (2.28) 

then escape from the trapping effect of the inversion. The ray with A equal to zero or 

0(s) = © is the limiting ray that has its turning point at infinity (the i = 2 case), while the 

• rays with i =3 correspond to negative values of A or Cos 0(s) > Cos© (0(s) < ©), and 

are the rays trapped by the ground. Rays with initial downward slopes can be divided 
in a similar manner but in all cases at least one reflection occurs before the ray 
^ escapes the trapping effect of the inversion, is the limiting case, or becomes trapped by 

the inversion. 

Figure 2.3 is an example of the rays calculated from the above equations for the 
case of an inversion. 


i 
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3.0 LAPSE CASE SOLUTION 


The solution of the problem posed by equations (1 .2) to (1 .4) in case of AT > 0, a 
lapse condition, was undertaken first. The general approach used in both the lapse 
and inversion cases was to first separated out a sinusoidal time dependence from the 
pressure, and then to Hankel transform the governing equations with respect to the 
horizontal distance from the source, r, to reduce the number of independent variables 
to one, the vertical height, z. This reduces the governing equation for the transformed 
independent variable to an ordinary differential equation for which an approximate 

solution can be obtained. This solution contains the Hankel transform variable, p, 
which replaced the horizontal distance in the transformed governing equation. The 
transformed solution must then be inverse transformed to return to physical space. 
Because of the complexity of the solution this inverse transform can not be carried out 
exactly and either an approximate inversion must be used or the inversion must be 
carried out numerically. Both approaches were used in the lapse case. In both of 
these approaches it is necessary to to interpret the Hankel transform variable as a 
complex variable and to continue the solutions off the real axis for the transform 
variable. This is not an intuitive process as the physical interpretation of the transform 
variable is lost off the real axis. This process will be discussed in detail below. 

With the solution obtained for complex values of the transform variable attention 
will be turned to the inversion of the transformed solution. The methods used are the 
classical saddle point approach and a Fast Fourier Transform (FFT) based numerical 
method. These methods are describe in detail elsewhere and will be only described 
briefly here. Finally the results of these approaches will be described. 


3.1 Transformation and approximate solution 


The time dependence in the governing equation and boundary conditions, 
equation (1 .3) and (1 .4) can be removed by assuming 


p(z, r, t) = e 10 * G(z, r) 


(3.1) 

The Hankel transform or two-dimension Fourier transform for an axisymmetric function 
can be defined [6] as 


G(z,P) = J G(z, r) r J 0 (|3r) dr 
o 


and the inverse transform as 

oo 

G(z, r) = J G(z, p) p J 0 (pr) dp 
0 


(3.2) 


(3.3) 

The use of transform methods in solving partial differential equations arises from the 
fact that an appropriate transform will convert a particular type of derivatives into an 

algebraic term expressed in terms of the transform variable ((5 in (3.2) and (3.3)) in 
place of the original physical independent variable (r in (3.2) and (3.3)). Thus the 
number of independent variables in the partial differential equation will be reduced by 
one and the transform variable acts only as a parameter in the transformed solution. In 


the case of Hankel transforms the radial dependence in the Laplacian operator 
expressed in cylindrical coordinates for a axisymmetric function is converted to an 
algebraic term, see [6] for more details. Applying (3.2) to (1 .3) leads to 




1 +az 


AT 

+ az + — 




(3.4) 

The term on the right hand side represents the source. The homogeneous form of this 
equation would have a solution with an oscillating behavior if the term in square 
brackets was positive and an exponential behavior if it was negative. Thus the 
transition occurs when the term is zero or 



Notice that if 


(3.5) 


B = -SL 
K a 


1 +as 


1 + as + 


AT 


Cos 0(s) 


(3.6) 

then equations (3.5) and (2.14) are very similar and we can interpret (3.5) as giving the 
value of p that causes the turning point (the location of the transition from oscillatory to 


exponential behavior as well as a horizontal ray) to be located at the height z. But (3.6) 
remains to be be interpreted. From (2.1 ) and (2.2), however, one finds that p is equal to 
(co/a J) times the cosine of the angle that a ray, which was initially at an angle 0(s), 
makes in the limit as height tends to infinity. Thus just as we have used 0(s) to identify 

a ray we can also use the Hankel transform variable {3. This discussion emphasizes 
the close relationship between the model being developed and the ray description. 
These parallels will be also be pointed out below. 

If height is nondimensionalized in (3.4) using the scale of the temperature 

gradient, a, a uniformly valid approximate solution to the resulting equation can be 

obtained, using the method presented by Nayfah [7], for large values of co/(a„a) as 


where 


G = —A—. h 1 (q(z, p)) + h 2 (q(z, p)) 

yj 9 2 (z. P) yj 9 z (z, p) 


(3.7) 


g 2 (z, P) = (1 +az + ~) 



i 
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1 


1 AT 


,«"( 44 ) 


2 T 


V'-(-e) : 


1 +0 


(3.8) 


and 




(3.9) 


(3.10) 


9 z (z. P) = 


3g(z, p) 

dz 


(3.11) 


The modified Hankel functions h 1 and h 2 are defined in [8] by 
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h .©'(f) V H ?’( I 4 2 ) 

3 

(3.12) 
and 

h 2 («“(|) 5 2 H ( f(|^) 

3 

(3.13) 

A and B are constants that have to be determined. This solution is rather complex and 
is expressed in terms of unfamiliar functions. Several important features of this solution 
need to be discussed to understand it. 

The functions h^) and h 2 (£) have a complicated behavior [8]. For real values of 
the argument both 1^ and h 2 yield a complex results which is oscillatory with an 
algebraic decay of the amplitude for increasing magnitude of the argument. The 
function h 1 can be shown to represent downward propagating waves (for e‘ ,<0 * as used 
here) and h 2 upward propagating waves. The oscillatory behavior also occurs for h, 

when the phase of the argument is equal to 2 jc/ 3 and for h 2 when the phase of the 

argument is -2rc/3. When the phase of the argument is tc/ 3, hj decays exponentially 

and h 2 grows exponentially. When the phase is -tc/3, h 2 decays exponentially and 1^ 
grows exponentially. 

The solution of the point source problem is closely related to the plane wave 


20 


problem discussed in [9]. In wave problems where transforms have been used the 
inverse transform can generally be interpreted as a superposition of plane waves over 
a range of angles. In the case of the Hankel transforms used here the limits of 

integration on p are from zero to infinity and p can be interpreted as oVa^ Cos 0^ where 

6„ is the angle between a ray and the horizontal in the limit of height tending to infinity 

as given in (3.6). Thus a value of p of zero corresponds to a ray which is propagating 

vertically upward at infinite height, a value of p of co/a^ corresponds to a horizontal ray 

at infinite height, and a value of P of ova^ {[1 + az ] / [ 1 + az + AT/T^ ]} 1/2 corresponds to 

a ray with imaginary slope at infinity such that its turning point (point of horizontal slope) 
is located at the height z. Based on this description the waves group themselves into 
several different forms. 

The first group, 0 <, p < p 0 = co/a» {1 / [ 1 + AT/T.. ]} 1/2 , are waves with their turning 
points at or below the ground surface and thus are actually reflected at the surface. 

These waves plus their reflections constitute the first group. The wave with p = p 0 
grazes the surface and is the limiting ray between the reflected and refracted rays. The 

second group, p 0 £ p £ p z = co/a.. {[1 + az]/[1 + az + AT/T.. ]} 1/2 , are waves with a 
turning point above the ground and below the observer height z. The waves in this 
group consist of those leaving the source in the range of angles described by (3.6) and 

their continuation after they have been refracted upward. The third group has p z £ p £ 

P 8 = co/a» {[1 + as ] / [ 1 + as + AT/T„ ]} 1/2 . These waves have their turning point above 
the observer and below the source. At the observers location these waves should yield 


a exponentially decaying solution. The waves in this group consist of those leaving the 
source in the appropriate range of angles and their continuation following refraction 
upward. These three groups of rays can all be seen in a ray diagram for a point source 
and all initially are propagating downward from the source. In addition there are waves 

propagating upward initially. These are in the range 0 < p < P s but differ from the first 

groups in that they do not originate as downward propagating waves that are reflected 
or refracted upward. Thus the "reflection" coefficient is missing from these waves. 

In addition to the above groups that can be seen in a wave diagram for a point 
source, there are several types that are necessary for the superposition given by (3.3) 

where p ranges from zero to infinity, but do not occur in a ray diagram for a point 

source. Group five consists has P s < p < co/a*,, these waves have the turning point 

above the source. In addition there are waves with co/a*, < p, these have no physical 

interpretation and correspond to complex angles at infinity. 

With these concepts let us proceed to the mathematical solution to the problem. 

The function g 3/2 given in (3.8) contains four branch points, two at p = ± p z and two at p 

= ± co/a... The negative branch points are not significant and will not be discussed. On 

the real p axis, for 0 £ p < p 2 , g 3/2 (z, p) is real. For p 2 < p < co/a*, g 3/2 (z,p) is positive and 

imaginary, and for p > co/a* g 3/2 (z,p) has a phase of -n at p = co/a*andtendstoa 

phase of -jt/2 as P tends toward infinity. This is shown in Figure 3.1 . The branch line 

for g(z,p) = ( g^z.p^can be chosen to be on the line where the phase of g 3/2 (z,p) is 

-jc. This line extends from the first branch point at p 2 to the second at co/a* and 


encloses a small region above the positive real axis, see Figure 3.2. The branches 
chosen for g(z,P) yield a phase of zero for 0 <, p < p z , tc/3 for p z < p < co/a and varying 
from -2rc/3 to -tc/3 for p > co/a^. 

Now by considering three cases the various forms of the solution can be 
obtained. These are shown in Figure 3.3. The first case is a wave with the turning 

point below the surface, 0 < p < p 0 . The second has the turning point below the source 

but above the surface, p 0 £ p < p s . In this case if the receiver is below the turning point 

then p > p z , if it is above the turning point then p < p z . The third case has P s < P ^ co/a., 
and the turning point is above the source. Again if the turning point is above the 

receiver then p > p z , if the receiver is above the turning point then p < p z . Note that 

these waves do not appear in a point source ray diagram but are needed to complete 
the solution. 

The solutions corresponding to the cases given above require the determination 
of the constants in (3.7). To do this a set of conditions are required. As result of the 
source terms in (3.4) the solutions separate into at least two forms, one for the region 
below the source and one above. The radiation condition, requiring outgoing waves in 
the limit as height tends to infinity, requires A to be zero above the turning point for 
z > s. At the source height, z = s, the solution must be continuous 

lim G(z,P) = lim G(z,p) 

z->$ + z-*s. 

(3.14) 

and must satisfy 



.. 9G. .. 5G Q 

hm [— — ] - lim [— ] = ^- 
z-»s. dz z-*s+ dz 2n 


(3.15) 

which is obtained by integrating (3.4) from z = s - e to z = s + e and taking the limit e->0. 
At a turning point continuity is required, 


lim G(z,p) = lim G(z,p) 

Z— >Zjp . 2— >Ztp + 

(3.16) 

At the ground surface, z = 0, the required condition is the normal impedance condition 
(1.4) which can now be expressed as 


c _ iZ dG 
(Dp dz 

o 


(3.17) 


Using these conditions, equation (3.4), and the physical situations presented in 
Figure 3.2 the following solutions can be obtained. For z > s 


G = K h 2 (r|(z,p)) [ IMq&P) + R 0 h 2 (Ti(s,P)) ] 


(3.18) 


forco/a eo >p z >p s >p 0 >p>0, 
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G = K h 2 (ri( 2 ,p)) [ MrKs.p) + R, h 2 (n(s,p)) ] 


f o r ca/a^ > P 2 > P 8 > P > P Q > 0, 


(3.19) 


G = K h 2 (t|(z,P)) e iw/3 R, [ h^iK^p) + R 0 h 2 ft(s,P)) ] 


for co/a. > p 2 > p > p s > p Q > o and 


(3.20) 


G = K [ h 1 (T|(z,P)) e i **+ R 0 h 2 (ii(z,p)) ] e** 
Ri[h 1 (Ti(s,p) + R 0 h 2 (Ti(s,p))] 


for co/a.. > p > p 2 > p s > p D > 0. For z < s 


(3.21) 


G = K h 2 (ri(s,p)) [ MTKZ.P) + Ro h 2 (n(z,p)) ] 


forca/a 00 >p s >p 2 >p 0 >p>i 


(3.22) 


G = K h 2 (r|(s,P)) [ h 1 (ri(z > p) + R, h 2 (n(z,p)) ] 


for co/a„ > P 8 > P 2 > P > P 0 > 0, 


(3.23) 
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<3 = K h 2 (rj(s,P)) e iw3 R, [ h^nCz.p) + R 0 h 2 (rKz,p)) ] 


for ro/a.. > P s > p > p z > p e > 0 and 


G = K [ h^rKs.p)) e^+ R 0 h 2 (r|(s,p)) ] 
Ri[h 1 (Ti(z,p) + R 0 h 2 (Ti(z,p))] 


for o/a« > p > p s > p, > p. > o. Where 


K = K(g (z,p), gjs.p)) 1 


12 i X 3 * v 9 2 ( s »P) 


R o = ' 


^^(TiCO.P)) + i>h 1 ’(r|(0,p)) 
x h 2 (Tj(0,P)) + i y h '( tj(0,P)) 


t=c a-'-Z 9 “ (0 ’ p) 
2 p„ a . 9 2 (°.P) 


(3.27) 


26 


and 


V = ' 


Po a . 


<^)*g 2 (o.P) 


(3.28) 



e 


.it 

*6 


.n 



(3.29) 

The solutions obtained above are for real values of p, but as discussed above p 

* must be interpreted as a complex variable. The boundaries between solutions off the 

real axis must be chosen as the branch lines used for calculating the function g(z,p), 
g(s,p) and g(0,p) from g 3/2 (z,p), etc. as were discussed above. On crossing these 

* branch lines it should be noted that the phase of g(z,P), g(s,p) and g(0,p) 

discontinuously jumps from -2 ti/ 3 to 2 jc/ 3 and the phase of g 2 (z,p), g 2 (s,p) and g 2 (0,p) 
increases by -2rc/3 (since it contains the root of g in the denominator). 

* Reference [8] presents some results for the modified Hankel functions that are 
useful for this type of behavior: 
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and 


(3.30) 


il* 

h^rie 3 ) = -h 2 (n) 

(3.31) 

Using these relations the solution as given by equation (3.18) to (3.25) can be rewritten 
and the regions of validity determined. For z > s these are 


G = K(g z (z,P), g z (s,P)) {h 2 (Ti(z,0)) [ Mn&P)) + RJgfO.P)) h 2 (n(s,p)) ]} 
in region A of p*space as given in Figure 3.4 


(3.32) 


G = K(g z (z,p), g z (s,p)) {h 2 (r|(z,p)) [ h^nCs.p)) + Ro(g(0,p) e'***) h 2 (ri(s,p)) ]} 

(3.33) 


in region B 


G = K(g z (z,p)e- i27t/3 , g z (s,P)) {h 2 (ii(z,P)) [ h^s.p ) &***) 
+ R o (g(0,p) e' ) h 2 (T|(s,p) e* 2 ** ) ]} 


(3.34) 


in region C 


G - K(g 2 (z,P)e- 1 2 * 3 . g z (s.P) e- ' ) {h 2 (t|(z,0) s' *** ) [ h, (ri(s.p) e> 2«3 ) 

+ R,(g(O.P) e'*«« ) h 2 (ii(s,p) 0 ^)]) 


in region D. For z < s 


(3.35) 


G = K(g z (z,P), g z (s,p)) h 2 (Ti(s,P)) [ h^nC z,P)) + R 0 (g(0,p)) h 2 (Ti(z,p)) ] 
in region E 


(3.36) 


G = K(g z (z,p), g z (s,p)) h 2 (n(s,p)) [ h 1 (n(z,p)) + R,, (g(0,p) e< ** ) h 2 (itfz,p)) ] 


in region F 


(3.37) 


G = K(g z (z,P). g 2 (s,P) e* ‘ ) h 2 (n(s,P)) [ h^TKz.p) e' ) 

+ Ro (9(0, P) e- ^ ) h 2 (ii(z,P) e i2lc/ 3)] 


in region G and 


(3.38) 


G = K(g z (z,p)e* * 2 n^ ( g z (s,P) e- * ** ) h 2 (n(s,p) 4 *** ) [ h 1 (n(z,p) ) 

+ R 0 (9(0.P) ©' 2,1/3 ) h 2 (n(z,p) e'2>t/3 ) ] 


(3.39) 


in region H. To obtain these results it is necessary to recognize that 


R 1 (g(o.P)) * R o (g(0,P) e 



) 


(3.40) 

From these results and the description of the behavior of the function g(z,p) (and 

therefore 7 \) at the branch lines it should be clear that the solution is continuous at the 

branch lines even with g(z,P) being discontinuous. This transformed solution must 
now be inverted. 


3.2 Lapse Results 


0 


Two methods were used to approximately invert the Hankel transform contour 
integration which lead to the saddle point method in the insonified region of physical 
space and a FFT based numerical method. Neither of these methods will be described 
in detail as the basic method is well known in both cases and the specific application 
has been described in detail elsewhere. 

Both of these approaches are based on the concept that although the inversion 
integral (3.3) is defined as being carried out along the real axis, the residue theorem of 
complex variables [10] allows the path of integration to be changed provided that there 
are no poles of the integrand between the original and modified paths. If poles exist 
then additional terms must be included with the integral along the modified path. In the 
case of an isothermal atmosphere the additional term due to the pole leads to the 
surface wave term. In the lapse case the only possible pole is the due to the 
denominator of the term multiplying the upward going wave (the reflection coefficient in 
the case where a reflection occurs) being equal to zero. These case has not been 

completely examined but in the limit of AT equal to zero it reproduces the surface wave 
term. Thus one clearly expects to see a similar behavior in the case of weak lapse 
condition. This surface wave like-behavior has not been investigated beyond the point 
described above and has not been included in the results given below. 

3.2.1 Contour Integration-Saddle Point Method 

Integrals of the form 
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+ 


l(5)=Je‘ 5W) dp 

c 

(3.41) 

can be approximately evaluated by the saddle point method if the path of integration C 
is such that the ends of the path do not significantly contribute to the integral, % is a 

large parameter, and f(p) has a point where the first derivative is zero. Complex 
variable theory indicates that a function can not have maximum in the region where it is 
analytic and a point where its first derivative is zero must be a saddle point [10]. At a 

saddle point if the path of integration follows the line of constant real part of f(P) then 

the imaginary part f(P) either increases or decreases at a maximum rate. If the path of 

integration follows the line of constant imaginary part of f(P) that passes through the 

saddle point then the real part of f(P) increases or decreases at a maximum rate. If we 

choose to follow a line of constant real part of f(p) through the saddle point in the 
direction such that the imaginary part increases at the maximum rate then the 
magnitude of the integrand decrease rapidly as we move away from the saddle point. If 

£ is large then the only significant part of the integral is near the saddle point and the 
integral can be approximated by using the first two non-zero terms in the Taylor series 

for f(p) yielding the well known results given in [1 1]. 

This method works well when a saddle point exists. However when one does not 
exist then an approximate integration can be carried out as described in detail by Ma 
[ 12 ]. 

To apply this method to the integrals given by (3.3) with (3.32) to (3.39) both the 
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Hankel functions and the modified Hankel functions must be replaced by their 
asymptotic expansions and the integral regrouped to extend from negative infinity to 
positive infinity. When this is done the result contains twenty terms since the integrand 

is different between in the various regions in {3-space and in each part of the G function 
contains two terms. Thus writing out only the terms of interest for z > s yields 


G(z,P) 


+ J lyz.P) e ' 1 P ‘ - a > d p 
0 


° V2 3/2 

xfir/ 7 RU' it?ra|9 (2.P) * 9 (s.to]-*/2} ft 
+ J K 7 (z,P) e dp 

Po 


+ ... 


jK 16 (z.p)e- UPr+Xt0 M) * fl M 


3/2 3/2 3/2 

‘ 2 g (0,f»]-n/2) 


dp 


Ps 


X. fk- /7 R\o' i,Pr+Xl9 (Z,P) + 

+ J K 17 (z,P)e 


3/2 3 12 

g (s,P) ] - it /6 } 


dP 


+ ... 

(3.42) 


< 
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and for z< s 




G(z,r) 


Po ^ 

+ 1 C 6 (z,r) e' ' ( P ' * M « «*•» • ”<*•» 1 ■ ■ ■ • ■ « > d p 
0 


. 3/2 3/2 

+ Jc 7 (z,P)e i|s '* M3 (s ' w -» w >*" ,a > d p 

Pc 


+ ... 


J 


3/2 3/2 3/2 

C 6 (z |3) e * i{pr + xts (z,p) + 9 < S ’W* 2 9 <o.P)]-*/ 2 ) d p 


* 3/2 3/2 

+ Jc i7 (z,p)e lni, * Ms w, * ( <s ' wl -" s) dp 

Po 


+ ... 


( 3 . 43 ) 


To find the saddle points the arguement of of exponential term must be differentiated 
with respect to (5 and set equal to zero. On differentiating one finds 
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9 r ^ 

3, a / 1 +as + y“ 

^(Xg 2 (z.P)) = -F(z,Cos' , (^\/ — P)) 

op (o V 1+as 

(3.44) 

where F(z,8(s)) is the function defined in (2.10) and used to describe the rays. 
Differentiating the arguments of the exponentials then yields equations (2.6) to (2.9), 
the equations defining the rays. Thus the saddle points associated with a particular 

point in physical space (in the insonified region) correspond to the values of p (or 0(s) 
where the two are related by (3.6)) that defines the two rays passing through the point 
in physical space. Just as the rays were interpreted as upward-going-direct waves, 
refracted waves, etc. the terms in (3.42) and (3.43) also have the same interpretations. 
Determination of the saddle point values then first requires determination of the types of 
waves present at a particular physical location and then solution of the appropriate two 
of equations (2.6) to (2.9). Once the location of the saddle point has been determined 
the classical results may be applied. A computer program for carrying out this 
procedure and the resulting equations to approximate the inversion integral have been 
given in detail by Cheng [13] and will not be repeated here. A typical result is shown in 
Figure 3.5. 

As the physical location of the receiver moves into the acoustic shadow real 

values of p or 0(s) cease to exist. Ma [12] has suggested an approximate approach 
which is also based on contour integration. In this approach the inversion integral 

between P 2 or p s and co/a.. is carried out numerically and the remainder of the integral 

extending from negative infinity and to positive infinity are carried out in a manner 
similar to the saddle point method. The integral carried out numerically physically 
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represents the contributions of the exponentially decaying disturbances due to waves 
with turning points above the receiver. Ma [12] also incorporated into his program 
Cheng's saddle point method for the insonified region with some changes. A typical 
result is shown in Figure 3.6. 

Both of these methods suffer from discontinuities as the receiver passes from a 
region where the receiver senses a direct and a reflected wave to one where a direct 
and a refracted wave would occur. This results not from the transformed solution which 
is continuous, but from the large argument approximation that must be made to obtain 
the saddle point form (3.41 ) and from the fact that the argument becomes zero at most 
of the boundaries. As a result of these intrinsic problems with the saddle point method 
a numerical approach was then applied. 

The saddle point method has the appeal of a physical interpretation of the 
mathematical steps and results. A purely numerical method loses that interpretation 
and the physical insight that comes from it. 

3.2.2 Numerical Integration Method 

The numerical approach used was developed by Richards and Attenborough [14] 
and was applied to the present case by Lloyd [1 5]. The method approximates the 
inversion integral by using a Fast Fourier Transform (FFT) algorithm. To obtain an 
integral suitable for the use of the FFT algorithm the Bessel function containing the 
horizontal distance dependence must be approximated by its asymptotic expansion. 
Three other modifications are then carried out. First, the integration path is modified to 
be above the real axis (Richards and Attenborough’s original approach was to 
integrate below the axis but they also assumed e-' 0 *.), this avoids the discontinuities at 
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the branch points but requires integration up the imaginary axis. Second, the integrand 
is modified to make the integral along the imaginary axis zero, but, as claimed by 
Richards and Attenborough, to not change the result. Finally an approximate term is 
added to account for the finite upper limit and to approximate the integral to infinity. 

This approach is described in detail by Lloyd [1 5] and Figure 3.7 is a typical 
result. Disadvantages of this approach are the large amount of computer time required 
and that a entire horizontal profile must be obtained at each receiver height. Thus to 
obtain a vertical profile many time consuming computer runs must be made and one 
point out of several thousand points is actually used from each run. This approach 
clearly does not contain the discontinuities present in the saddle point method. Figure 
3.8 compares the saddle point method and the purely numerical method. The 
agreement is excellent in the insonified region with the exception of the region very 
near the shadow boundary. The agreement is good in the initial sound level decrease 
as the shadow boundary is crossed but the saturation region deep in the shadow is not 
the same for the two methods. 

The numerical method often results in oscillations in the sound level at large 
distances from the source, this appears to be an artifact of the numerical inversion 
method and is dependent on the parameters of the inversion scheme. Also as very 
large distances are approached the calculated sound level often increases this is 
clearly due to the numerical inversion method. These points are further discussed by 
Lloyd [15] and a listing of Lloyd’s program is given in Appendix A. 


37 


4.0 INVERSION SOLUTION 


The inversion solution, AT < 0, follows very closely along the lines of the lapse 
solution. The governing equation and boundary conditions, (1.2) to (1.4), are identical 
in the two cases and the approach using Hankel transforms is also the same. Again the 

solution requires analytic continuation off the real p axis. This has proved to be difficult 
and lead to several errors initially. (The solution given in the Sixth Semiannual Report 
[18] are incorrect.) 

Although the solution can be discussed in terms of rays and the closely related 
saddle point method, this approach has not been used to approximately evaluate the 
solution in the inversion case. In the lapse case only two rays, at most, pass through a 
given point. In the inversion case, at large distances from the source, many rays may 
pass through a given point due to the "trapping" effect of the inversion. Since the 
saddle point method requires all of these rays and their corresponding saddle points to 
be located, and this is the most difficult part of the method, the approach becomes 
impractical. Thus only the purely numerical method of inversion has been used. 

4.1 Transformation and approximate solution 

The time dependence is removed from (1.3) as in (3.1) and the resulting equation 
Hankel transformed using (3.2) to obtain (3.4). Again the location (in terms of the 

transform variable P) of the transition from oscillating to exponential behavior is given 

by (3.5). However, since AT/T.. < 0 the transition is at a value of p greater then co/a... 
Using (3.6) and comparing results to those of Section 2.2 one can interpret the solution 


in the range 0 £ J3 < co/a.. as representing the rays that will escape the trapping effect of 

the inversion. These rays either go directly from the source to infinite heights or go 
from the source to the ground where a reflection occurs and then go to infinite heights. 

In the range co/a.. £ p £ p s the solution represents rays that are trapped by the 

inversion. The transition given by (3.5) applies to rays in this range. Beyond this range 
the rays do not occur in a ray diagram. 

The solution of (3.4) can again be approximated by (3.7) and (3.9) through (3.1 1 ). 
Equation (3.8) must modified by a negative sign on the right hand side yielding 


3 

2, Q . . 4 AT . 

g (z,(3) = - ( 1 + az +— ) 




1 +<E 
1 -<D 


(4.1) 


This change is necessary since the region of oscillatory solution is below the turning 
point in the inversion case while it was above it in the lapse case (see Nayfeh [7]). In 

the inversion case AT/T^ < 0 and thus <X> > 1 for 0 < p < co/a It is convenient to note 
that we may rewrite the logarithm term in this case as 


to clearly indicate the choice of the branch of logarithm, ln(-1) = +i it, the opposite of that 
in the lapse case. Thus, in the range 0 <, p < co/a« g 3/2 (z,p) ranges from slightly below 
the negative real axis to infinity along the negative imaginary axis. For values of P 
between co/a« and P 2 ranges along the positive real axis from infinity to zero. For real p 

greater than P 2 values of g 3/2 (z,P) are along the positive imaginary axis, see Figure 4.1. 
As in the lapse case the boundary between these regions are branch points of the 
function g 3/2 (z,p) with the branch lines extending downward from the branch points for 
Re(p) > 0 in p-space. 

The argument of the Hankel functions involves g(z,p) = (g 3/2 (z,P)) 2/3 and again the 
branches must be chosen with care. For 0 £ p < co/a^ g(z,p) is chosen such that its 
phase ranges from slightly greater then zero (or 2%) at p = 0 to it/3 as the branch point 

at p = (o/a.. is approached from values of p less then oVa.. In this region the two 
modified Hankel functions have an oscillatory and exponential growth or decay 
behavior with with one (h-,) representing upward traveling and decaying waves and 

the other (h 2 ) downward traveling, growing waves. In the range co/a*. < p < p 2 and 

g 3/2 (z,p) is real and positive. The function g(z,p) is chosen to be on the line with phase 


-2jt/3. In this case the modified Hankel function h 1 represents both upward and 
downward traveling waves, a standing wave-like phenomena, and h 2 a downward 
traveling wave. 

The branch line for g(z,p) = (g 3/2 (z,P)) 2/3 needs to be defined to extend these 

solutions off of the real axis. Lines of constant phase of g 3/2 (z,P) run from the first 
branch point to the second in a manner similar to the lapse case, but with the order of 
the branch points reversed. Figure 4.2 shows the behavior of these lines of constant 
phase. Again a line of constant phase is a convenient branch line. 

If the line where the phase of g 3/2 (z,{3) equals -k/2 is chosen as the branch line for 
g(z,p) = (g 3/2 (z,p)) 2/3 then the phase of g(z,p) can be made to agree with the desired 
values on the real axis as described above. In addition for real p and p z < p, g(z,p) 

has a constant phase of it/ 3 with h 1 (Tt(z,p)) having a decaying exponential behavior for 
increasing z and representing the contribution of waves with a turning point below the 
receiver’s height to the total pressure field. 

Using the conditions given in (3.14) to (3.1 7) and the physical descriptions of the 
type of waves that occur in each situation the following solutions can be obtained. For 
z>s 

G = K [ iyrKs.P)) + R 0 h 2 (n(s,p)) ] h 2 ft(z,P)) 

(4.3) 


which is identical to (3.18) for p 0 > p s > p 2 > to/a 00 > p> 0, 


G«-K[ h,(n(»,p)) / R„ + h 2 (n(s,p» ] h,(Ti( 2 .P) 


tafP.>p t > Pj>pxo/a.>0, 

.2r 

G = ■ K e 3 [ h,(n(s,p)) / R 0 ♦ h 2 (n(s,p)) ] h 2 (n(z,p) 

fo r P 0 >P s >P> P z >co/a 00 >Oand 

G = K [ h,(7l(S,P)) + R 2 hj(n(8,P)) 1 h 2 (r,(2.p)) 


forp 0 >p> p 8 >p 2 >co/a M >0. For s > 2 

G = K [ h,(n( 2 ,p)) + R 0 h 2 (t|(z,P)) ] h 2 (n(s,p)) 


the same as (3.22) for p Q > p z > P s > oVa_ > p > 0, 


G = - K [ h,(ii(z,p)) / R 0 + h 2 (n(z.p» ) h,(n(s,p) 


for P Q > P 2 > P, > P > oVa_ > 0, 




.2k 

*T 

G = - K e [ h 1 (Ti(z,P)) / R o + h 2 (t|(z,P)) ] h 2 (ri(s,P) 

for P 0 > P 2 > P > P s >Gl/a ~ > 0 and 


(4.9) 


G = K [ h 1 fo(z.p)) + R 2 h 2 (ri(z.p)) ] h 2 (Tj(s,p)> 


(4.10) 

P c > P> P 2 >P S > Here K and R Q , and R 1 are defined by (3.25) and (3.26) 

and 


R 


2 



( 


R 


o 



R 


o 


) 


(4.11) 


As discussed above these solutions must be continued off the real axis. As in the 
lapse case the boundaries are chosen as the branch lines for calculating g(z,p), g(s,p) 
and g(0,p) from g 3/2 (z,p), etc. On crossing these branch lines the phase of g(z,p), g(s,p) 

and g(0,p> jumps discontinuously from ji/ 3 to - n and g z (z,p), g 2 (s,p) and g z (0,p) jumps 

discontinuously by 2 tc/ 3. Using (3.30) and (3.31) solutions (4.3) through (4.10) can be 
continued off the axis as 


l 
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G = K(g z (z.p), g z (s,p)) [ h 1 (n(s,p)) + R o (g(0,p)) h 2 (n(s,P)) ] h 2 ft(z,p)) 


in region A of Figure 4.3 


j 2rt . 2n jiZL 

G = K(g 2 (z,p) e 3 , g 2 (s,p) e 3 ) [ h,(n(s,p) e 3 ) 

I . An . 4n 

+ R o (g(0,P) e 3 ) h (Ti(s,p) e 3 )] h 9 (n(z,p) e 3 ) 


in region B, 


I 2ft . 4 ji 

G = K(g z (z,P), g z (s,P) e 3 ) [ h 1 (r t (s,p) e 3 ) 

. 4 n . 4 n 

T 'T 

+ R o (9(0,P)e )h 2 (rKs,P)e )] h 2 (n(z,p)) 


in region C and 


. 4n 


(4.12) 


(4.13) 


(4.14) 


G - K(g z (z.p), g 2 (s,P)) [ h,(g(s,p)) + R o (g(0,p) e ) h 2 (n(s,p)) ] h 2 (n(z.p» 


(4.15) 


in region D. For z < s 


G = K(g 2 (z,P), g z (s,p)) [ h^itfz.p)) + R o (g(0,p)) h 2 (ti(z,p)) ] h 2 (r|(s,p)) 


(4.16) 
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in region E, 


- I- 


. 2 * 


• i- 


4* 


G = K(g z (z,p) e , g 2 (s,p) e ) [ h 1 (n(z.p) e ) 


. 4* 


. An 


[ , An 

* R 0 (9(°.P) e 3 )h 2 (n(2,|3)e 3 )] h a (T»(s.p) e 3 ) 


(4.17) 


in region F, 


i — | 4 k 

G = K(g 2 (z.P) e 3 , g z (s,p) ) [ h,<Ti(z,p) e 3 ) 

i~ iii 

+ R o (g(0,P) e 3 ) h.(n(s,P))] h.(r,(z,p) e 3 ) 


in region Q and 


(4.18) 


. An 

G = K(g z (z,p), g 2 (s,p)) [ h 1 (ri(z I p)) + R o (g(0,p) e 3 ) h 2 (T](z,p)) ] h 2 (Ti(s,P)) 

(4.19) 

in region H. 

From these results and the description of the behavior of the function g(z,p) (and 

therefore t|) at the branch lines it should be clear that the solution is continuous at the 

branch lines even with g(z,P) being discontinuous. The transformed solution must now 
be inverted using the numerical method developed by Richards and Attenborough [14]. 


4.2 Inversion Case Results 


As was discussed above the numerical method originally developed by Richards 
and Attenborough [14] was used to invert the Hankel transformed solution in the case 
of the inversion. This was due to the fact that many rays pass through a given point in 
physical space in the case of an inversion and the number of saddle points which exist 
equals the number of rays passing though that point. Since finding the saddle points is 
the most difficult and time consuming part of that method the approach appeared 
impractical in this case. The numerical method used is identical to that of the lapse 
case as described by Lloyd [15]. 

Only a limited number of cases have been mn to date using the solution 
described in Section 4.1 and the numerical inversion technique. Figure 4.4 shows a 
typical case. The results generally show an interference pattern with 6 dB/doubling of 
distance decay out to distances of the order of thirty meters and a more complicated 
behavior beyond that distance but with no significant change in the rate of decay. This 
latter result is somewhat unexpected from qualitative arguments. Experimental data for 
propagation under inversion conditions is quite limited, with the data presented by 
Sutherland and Brown [16] being the major set. However, this set contains only seven 
measurements at a fixed height over a 675 meter distance. No direct comparisons 
have been made but the data also shows what appears to be a 6 dB/doubling of 
distance decay with some interference minima. Thus at least qualitatively the 
agreement appears good. 

Appendix B contains a listing of the program for the Inversion case. 


5.0 CONCLUSIONS 
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Approximate solutions of the Hankel transformed acoustic wave equation with a 
particular, realistic and well-developed vertical sound speed (or temperature) profile 
have been obtained for both the lapse and inversion cases. These solution are quite 
complex and exact inversion of the transformed solution does not appear possible. 

Both approximate inversion using contour integration and the saddle point approach 
and numerical inversion have been used to obtain the physical solution in the case of 
lapse conditions. Only numerical methods have been used with inversion condition. 

The lapse case shows the expected behavior: an interference pattern with a 6 
dB/doubling of distance decay within the shadow region; a rapid decrease in sound 
level in the vicinity of the geometric shadow boundary; and approximately a 6 
dB/doubling of distance decay well within the shadow region. Similar behaviors occur 
for both inversion methods but the contour integration - saddle point method yields and 
larger decrease in the sound level on passing into the shadow than the numerical 
integration technique. The origin of this difference has not been determined. The 
contour integration - saddle point method results appears to agree with the empirical 
model of Weiner and Keast [17] better then the results of the numerical inversion 
technique. Since the techniques are applied to the same approximate solution of the 
transformed acoustic wave equation the difference must result from the inversion 
techniques. The numerical technique also produces a weak interference-like behavior 
far into the shadow region. This appears to be artifact of the numerical method as is the 
increase in sound level that frequently occurs as the maximum distance for the 
inversion technique is approached. 

Agreement between the results and experimental data is fair within the shadow 
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boundary. The level Is predicted well but the location of interference maxima and 
minima are not accurately predicted. This may be due to the poor fit of the temperature 
profile to the measured profile. No data appears to be available that both gives a 
temperature profiles and sound levels in the shadow region. 

The inversion case shows an interference pattern with a 6 dB/doubling of 
distance decay out to distances of the order of thirty meters for realistic temperature 
profiles. Beyond this distance the decay rate appears to remain nearly the same but 
the structure of minima and maxima becomes irregular. This tends to agree with a 
simple geometric argument since "trapped" rays start to reappear in a ray diagram at 
such distances. Little data is available for comparison in the inversion case. 



6.0 REFERENCES 


1 . Best, A. C., "Transfer of Heat and Momentum in the Lowest Layers of the 
Atmosphere," Meterological Office Geophysical Memorandum No. 65, 1 935. 

2. Geiger, R., The Climate Near the Ground . Harvard Univ., Cambridge, 1965. 

3. Reynolds, D. D., Engineering Principles of Acoustics . Allyn and Bacon, 
Boston, 1980. 

4. Butterworth, J. R., "An Investigation of the Sound Pressure Levels in a 
Acoustic Shadow," M.S. Thesis, University of Utah, 1979. 

5. Van Moorhem, W. K., Ma, Y., and Brown, J. M., "Ray Paths near the Ground 
in a Realistic Thermally Stratified Atmosphere," J. Acoustical Society of 
America, 80, 650, 1986. 

6. Bracewell, R. N., The Fourier Transform and its Applications . 

McGraw-Hill, New York, 1986. 

7. Nayfeh, A., Perturbation Methods . Wiley, New York, 1 973. 

8. Tables of the Modified Hankel Function of order One-Third and of their 
Derivatives . Harvard Univ., Cambridge, 1945. 

9. Van Moorhem, W. K. and Landheim, G. K., "The Propagation of Plane Waves 
in a Thermally Stratified Atmosphere," J. Acoustical Society of America, 

76, 867, 1984. 

10. Churchill, R. V., Complex Variables and Applications . McGraw-Hill, New 
York, 1960. 

1 1 . Mathews, J. and Walker, R. L., Mathematical Methods of Physics . 


Benjamin Cummings, Menlo Park, CA, 1964. 

1 2. Ma, Y. "Acoustic Wave Propagation Along a Finite Impedance Ground 
Under Temperature Lapse Conditions," M.S. Thesis, University of Utah, 

1984. 

13. Cheng, A. S. C., "Sound Propagation from a Point Source in a Thermally 
Stratified Medium," M.S. Thesis, University of Utah, 1985. 

14. Richards, T. L. and Attenborough K., "Accurate FFT -Based Hankel 
Transforms for Prediction of Outdoor Sound Propagation," J. Sound and 
Vibration, 109, 157, 1986. 

1 5. Lloyd, S. R., "Nonisothermal Acoustic Propagation Analysis using Fast 
Fourier Transforms," M.S. Thesis, University of Utah, 1985. 

16. Sutherland, L. C., and Brown, R., "Static Tests of Excess Ground 
Attenuation at Wallops Flight Center," NASA Contractors Report 3435, 

1981. 

17. Weiner, R. M., and Keast, D. N., "Experimental Study of The Propagation 
of Sound over Ground," J. Acoustical Society of America, 31, 724, 1959. 

18. Van Moorhem, W. K., "Sixth Semi-Annual Report to the National 
Aeronautics and Space Administration on Grant NAG-1 -283 Acoustic 
Propagation in a Thermally Stratified Atmosphere," UTEC ME 86-034, 
Mechanical and Industrial Engineering Department, University of Utah, Salt 
Lake City, 1 986. 


7.0 LIST OF SYMBOLS 


English 



Sound speed. 

Function defined by (2.4) or constant in (3.7). 

Function defined by (2.5) or constant in (3.7). 

Function defined by (2.6). 

Constants. 

Function defined by (2.7). 

Function defined by (2.1 1 ). 

Arbitary function 
Function defined by (2.10). 

Function defined by (2.25). 

Function defined by (2.26). 

Function defined by (2.27). 

Function defined by (3.8). 

Hankel transform of G. 

Acoustic pressure with time dependence seperated out, see 
(3.1) 

Modified one-third order Hankel function of the first kind, see 
(3.12). 

Modified one-third order Hankel function of the second kind, see 
(3.12). 

V-1 

Intergal defined by (3.41). 

Zero order Bessel Function. 

Function defined by (3.25). 

Constants 

Acoustic pressure. 

Constant determining the strength of a point acoustic source. 
Horizontal distance from the source. 


r Horizontal distance from the source. 

R 0 Function defined by (3.26). 

R, Function defined by (3.29). 

Rg Function defined by (4.1 1). 

s Height of the point source above the ground, 

t Time. 

T Temperature, 

z Height above the ground surface. 

Z Acoustic impedance of the ground surface. 


Greek 

a Scaie factor for temperature, see (1 .2). 

p Hankel transform variable replacing r, see (3.2). 

Y Function defined by (2.1 7). 

5 Delta function. 

t\ Function defined by (3.10). 

0 Angle an acoustic ray make relative to horizontal. 

© Limiting ray angle, see (2.28). 

X aj(co a). 

£ Arbitrary arguement 

p Density of the air. 

x Function defined by (3.27). 

4 Function defined by (2.12). 

4> Function defined by (3.9). 

H* Function defined by (3.29). 

co Circular frequency. 
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Combinations 


AT Change in Temperature between the ground surface and 

far above the ground. 

Subscripts other than given above 


i 

1 ,2 or 3. 

0 

Evaluated 

tp 

Evaluated 

s 

Evaluated 

z 

Derivative 

oo 

Evaluated 


at the ground, or a reference value, 
at a ray turning point, 
at the source height s. 

with respect to height (g z or g z2 ) or evaluated at the height z. 
at infinite height. 
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Figure 1.1. Temperature as a function of height above the ground for different times of 
the day as determined by Best [1]. 
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Figure 1 .2. The present model of temperature as a function of height and a set of 
observations by Butterworth [41. 
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ELEVATION 



Figure 2.2. Acoustic rays for a lapse case with a = 1 .75 m* 1 and AT/T^ = 0.03 with a 
source at a height of 2 m. 


57 



elevation (m> 



Figure 2.3. Acoustic rays for an inversion case with a = 1 .75 m' 1 and AT/T^ = -0.03 
with a source at a height of 2 m. 
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V 1+oa+ t: 



Branch lines 
for g 3/2 





Figure 3.2. Sketch of the location of the branch line used for calculating g (z,|3) 
( g 3/2 (z,p) )2/3 in the lapse case. 




Figure 3.3. Sketch of the three types of physically occuring rays in the lapse case. 
Type 1 rays have their turning point below the ground surface. The turning point for 
type 2 rays is below the source and above the ground surface. Type 3 waves have a 
• turning point above the source, this type of ray does not appear in a point source ray 

diagram but occurs in the superposition making up the inverse transform. 
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Figure 3.4. Sketchs of the regions in complex (5-space where the various forms of the 
solution are valid for the lapse case. Part a) is for points above the source, part b) is for 

points below the source. The lines are branch lines for g (0,P), g (z,P) and g (s,P). 
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SPL-dB 



Figure 3.5. Typical results for a lapse case using the saddle point method only, from 
Cheng [13], as compared to the Weiner and Keast empirical model [17]. The solution 
extends only to the shadow boundary at about 68 meters. 
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Results of Ma{12) 


Results of Weiner 
and Keast[17] 


Excess Attenuation 

s * 3 meters \ 

z « 1 meter \ » 

© m 1 0,000 rad/sec (1 590 Hz) ^ ^ 

a = 3.5 (meters)* 1 • / ; 

AT/T„« .001 V 


r - meters 


1000 


Figure 3.6. Typical results for a lapse case using the combined saddle point-contour 
integration method. From Ma [121. 



Sound pressure level (dB) 



Figure 3.7. Typical results for a laps© case using the nurnarical inversion technique. 
From Lloyd [15]. 
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z » 1 .5 meters 

co = 10,000 rad/sec (1590 Hz) 
a = 2.5 (meters)* 1 
AT/T 0 .= .025 
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and Keast[17] 

Analytical 

approximation^] 

Numerical Result[l5] 
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Figure 3.8. Comparison of the results of the saddle point-contour integration method 
and the numerical technique for a lapse case. 


•T* 





Figure 4.1 . Path followed in complex g3/2 . space as the rea( part of p varies from 2ero 
to infinity and the imaginary part of p is constant for the inversion case. 
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Im(p) 


• 9 




Figure 4.2. Sketch of the location of the branch line for calculating g (z,{3) = 
( g 3/2 (z,p) )2/3 in the inversion case. 
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Figure 4.4. Typical results for an inversion case using the numerical inversion 
technique. 
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APPENDIX A - The Lapse Case Program 


The material below is from Lloyd [1 5] and both describes the program used to 
calculate the transformed solution and to carry out the inversion and presents a listing 
of that program. The program is named FFTPRESS and was written to be used on a 
VAX 750. Descriptions are given both of the subroutines comprising the program, of 
the input variables and a set of "helpful hints" are given that may be useful in running 
the program. The intent of this section is not to describe in detail the operation of the 
program but to allow a somewhat experienced Fortran user to run the code as it was 
created. 

SUBROUTINES 

Input: Subroutine to input the necessary parameters to the main 
program. The following sentences summarize each of the input variables 
in the order they are requested. Tint is the temperature at infinite height, 
normally 300 K. Tinf is used to calculate the speed of sound a. Dtot is the 
temperature change from infinity to the ground normalized by the 
temperature at infinity. Dtot is normally 0.025. Alpha is the term used in 
the temperature profile defining the altitude at which the temperature 
gradient becomes effective. Alpha is normally 2.5 (meters)' 1 . Splref is 
the reference sound pressure level used in the calculations of sound 
pressure level in dB. Omega is the frequency of the sound source in 
rad/sec. Resistance is the flow resistance used in the Chessel model and 
is normally 300 cgs units. Zr is the height of the receiver in meters. Zs is 
the height of the source in meters. Alp is the amplitude of the imaginary 
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component in transformed space. Alp greater than zero is integration 
above the real axis. The product of Alp and the step size AK should not 
exceed 0.01 and may be much smaller. Alp is replace by |3 in the thesis. Ns 
is used in an analytical function that is subtracted from the sampled 
solution to null the effects of off axis integration. Ns is normally 3. If Ns 
is equal to zero then no subtraction occurs. Ns is replaced by 8 in the 
thesis. Me nonzero signals the inverse Hankel transofrm routine that the 
terms representing the integral extended to infinity are to be included in 
the inversion. Me is also the number of terms to be used and is normally 5. 
Me is replaced by M in the thesis. N1 is the number of points to be used. 

N1 depends on the maximum horizontal distance desired. N1 equal to 4096 
points is a common value. Np and N1 are used interchangeably. N1 must be 
equal to an integer power of two. Delbeta or delK are the step size in 
complex K space. In the program Delbeta and DelK are used 
interchangeably. Delbeta also depends on the maximum horizontal distance 
desired and also on the maximum Beta allowed. This maximum Beta is very 
near to omega divided by the speed of sound. Beta and K are used 
interchangeably in the program and thesis. 

Region: Subroutine used to determine which of 8 diffferent forms of 
the general solution are to be used. The selection depends on how the 
waves are interferring at that particular value of K. Region calls to 



subroutine g32all to identify if the complex part of g 3 ^ 2 has changed sign. 
The sign change indicates a different set of rays are combining to form the 
solution. With each set of rays a different form of the solution is required. 

Hall: Subroutine to calculate the Hankel functions H-j and H 2 and their 

derivatives in terms of the Airy functions, Al and Bl. Hall calls to Cgbair 
to get the Airy functions needed. 

Hall2: Subroutine to calculate only the Hankel functions. 

Cgbair: Subroutine to calculate the Airy functions. Cgbair uses either 
an asymptotic or a small argument approximation of Al and Bl depending on 
the value ofcomplex K. 

Gzalll: Subroutine to calculate the derivatives of the g 3/2 function. 
These values are used to compute the reflection coefficients. 

Gall: Calculates the g function needed to calculate the g 3/2 function. 

Dafb2: Subroutine modified from Attenborough and Richards to 
calculate the inverse Hankel transform using fast Fourier transforms. 

Dafb2 calls to subroutine Zeta to calculate the terms representing the 


extension of the integral to infinity and to subroutine Fork used to perform 
the actual fast Fourier transform. 


Zeta: Subroutine to compute the value of the integral to infinity. 

Fork: An efficient fast Fourier transform taken from Attenborough 
and Richards program. 

VARIABLES 

Beta: Used interchangeably with K, both are the complex argument of 
the transformed solution. 

Tau: The term i used in the reflection coefficients developed by Van 
Moorhem. 

Sci: The term y used in the reflection coefficient developed by Van 
Moorhem. 


Rq: The actual reflection coefficient. 


R-j : The modified reflection coefficient representing refraction. 




En: The complex argument of the Hankel function. En is a function of 
K (or Beta) and the height of the source or receiver, whichever applies. 

Rlemda: The ratio of omega to speed times alpha. The reciprocal of 
the wave number. 

Zimped: The complex impedance of the ground normalized by the speed 
of sound and the density. 

Z2 and Z3: Heights of the receiver and source, respectively. 

Gbar (K) or Gbar(Beta): The sampled function to be inverted. 

Gbar(r): The inverted solution. The real space answer. 

G: The sound pressure level result. 

Rad: The horizontal distance of the present (Gbar (r) 

Rad2: The logarithm of Rad. 
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HELPFUL NOTES 

1 . All units used here are examples only. The only requirement in 
operation is to use a consistent set of units. 

2. The output of sound pressure level and horizontal distance along 
with an echo of input parameters is written to file FOR044.dat in the 
present diretory. 

3. To plot the result at the University of Utah Mechanical Engineering 
Vax system type RUN PLOTPROMPT and the rest is interactive. 

4. Basic instructions for use on the University of Utah Mechanical 
Engineering Vax system are: 

a. Log on using normal sequency of user name and password. 

b. Type @Q to link all necessary files together. Instead of combining 
all files into a large file several small trackable files are used for ease of 
editing. 

c. Type RUN FFTPRESS to begin execution. 

d. Input the variables as requested by subroutine Input. 

e. At completion FFTPRESS will display FINALLY FINISHED. To plot the 
results type RUN PLOTPROMPT. This is a standard plotting program that 
uses the system subroutine Mgraph. The data file name is FOR044.DAT. 

The data file contains 2 columns. The fist column of the data is the 
logarithm of the horizontal distance. The second column of the data is the 
sound pressure level. 15 lines are used at the beginning of FOR044.DAT to 


echo the input parameters, therefore tell PLOTPRMPT to begin accepting 
data at line sixteen of file FOR044.DAT. The plot will display on graphics 
terminals only. Mgraph asks if a hard copy is desired when crt plotting is 
finished. PLOTPROMT has autoscaling capability that can be turned on or 
off and offers many other self instructing options. Mgraph creates files 
named HPPLOT.HPL, however, it is recommended to change the name as soon 
as possible to avoid deletion of previous plots. If a hard copy plot is 
desired after exiting PLOTPROMT type PLOT then the file name, 
f. Log off with command LO. 
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ooooooooooooooooooooooooooooooooooooooooooooooooooo 


PROGRAM FFT PRESS 


MAIN ROUTINE TO DEVELOPE THE GBAR(b,Z) TO BE INVERTED * 
THE MAIN WILL INPUT THE NECESSARY CONSTANTS, CALCULATE BETA, * 
SELECT FORM OF SOLUTION CALL THE SUBROUTINES AND COMPUTE RO,TAU, * 
SCI, g, THE HANKEL FUNCTIONS AND FINALLY CALCULATE GBAR(b,Z). * 

HHHHHHHHHHt 


SUBROUTINES: 

G32ALL: 

REGION: 

GZALL: 

GALL: 

HALL: 

HALL2: 

DAFB2: 


ZETA: 

FORK: 

INPUTS: 

TINF: 

DTOT: 


ALPHA 4 : 

SPEED: 

OMEGA: 

Z: 

S: 

SPLREF: 

RESISTANCE: 

ALP: 

NS: 

ME: 

N1 : 

DELBETA: 


FINDS THE FUNCTION g*3/2(b,z) USED TO DETERMINE 
WHICH REGION OF SPACE PRESENT BETA IS IN. 

GIVEN g*3/2 DETERMINES WHICH REGION OF 
SPACE PRESENT BETA 1$ IN. 

FINDS gz(t>,Z)=dg/dZ AND THE OTHER DERIVATIVES 
USED TO FIND K,TAU,SCI. ALSO RETURNS ETA(b,Z) 
THAT IS -USD AS THE ARGUMENT FOR HANKEL 
FUNCTIONS.' GZALL CALLS TO GALL TO GET g. 

CALCULATES g FUNCTION. 

CALCULATES THE HANKEL FUNCTIONS FROM AIRY 
FUNCTIONS. CALLS CGBAIR TO GET AIRY FUNCTIONS 

CALCULATES ONLY THE HANKEL FUNCTIONS NOT THE 
DERIVATIVES AS HALL DOES. CAMS CGBAIR ALSO. 

GIVEN GBAR(b,Z) USES METHOD DEVELOPED BY 
RICHARDS AND ATTENBOROUGH TO PERFORM THE 
HANKEL INVERSION. USES SUBROUTINES ZETA, FORK. 
DAFB2 MAKES SEVERAL CORRECTIONS TO A GENERAL 
FFT. THE STANDARD CODE IS TAKEN FROM RICHARDS 
AND ATTENBOROUGH PROGRAM. 

FORMS THE SUM OF ME TERMS WHICH APPROXIMATES 
GBAR(b,Z) TO INFINITY IN THE BETA SPACE. 

A VERY FAST FFT USED TO PERFORM ACTUAL 
INVERSION OF GBAR(b,Z) FROM THE BETA SPACE. 

THE TEMPERATURE AT VERY LARGE Z 

THE DELTA_T/T PARAMETER REPRESENTING THE 
TEMPERATURE GRADIENT. 

PARAMETER USED IN DEFINITION OF TEMPERATURE 
GRADIENT 

SPEED OF SOUND AT T(INF). 

FREQUENCY OF SOUND IN RADIANS/SEC. 

FIXED DISTANCE TO THE OBSERVER Z2 IN PROGRAM. 

FIXED DISTANCE TO THE SOURCE Z3 IN PROGRAM. 

REFERENCE SOUND PRESSURE LEVEL USED TO COMPUTE 

THE SOURCE STRENGH Q. 

GROUND RESISTANCE IN THE CHESSELL MODEL. 

THE TERM USED TO INTEGRATE OFF REAL BETA AXIS. 

PARAMETER IN THE ANALYTICAL FUNCTION IN DAFB2. 

PARAMETER TO PRODUCE SUM TO INFINITY IN DAFB2. 

SIZE OF ARRAY TO BE INVERTED. 

STEP SIZE USED FOR BETA AISO DELK IN DAFB2. 


NOTE: 


K AND BETA AND DELK AND DELBETA ARE USED INTERCHANGEABLY 


78 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


VARIABLES: 

BETA: INDEPENDENT VARIABLE IN TRANSFORMED SPACE. 

TAU: TAU DEFINED IN PAPER BY VAN MOORHEN. 

SCI: SCI DEFINED IN PAPER BY VAN MOORHEM. 

RO: RO CONSTANT DEFINED IN VAN MOORHEN. 

R1 : R1 CONSTANT DEFINED IN VAN MOORHEN. 

IN: HANKEL FUNCTION ARGUMENT. 

RLEMDA: OMEGA/ (ALFHA*SPEED) . 

ERO: BRANCH CUTS IN THE BETA SPACE ASSOCIATED WITH 

BRZ: THE SQUARE ROOT AND 3/2 POWER FUNCTIONS 

BRS: IN g AND g THREEHALF. 

BRW: BRANCH CUT AT OMEGA/SPEED. 

Z IMPED: GROUND IMPEDENCE NORMALIZED BY DENSITY AND SPEED. 

GZ: dg/dZ FROM SZALL1 . 

GZZ: d2g/dZ2 FROM GZALL1 '. 

OTHER DERIVATIVES £eR THIS NOTATION 
Z1 : REFERENCE DISTANCE 0.0. 

Z2: Z AS ABOVE. 

Z3: S AS ABOVE. 

I KXXXXXXX X XXEXX E XXXXXKXXXXXX X XKXXXXXXE X XXXXKXXXXX g X X -E* E-EHHt-E E E-E-E-E- E -E- K -E 

IMPLICIT DOUBLE PRECISI0N(A-H,0-Z) 

INTEGER IREGION 
COMMON /INTBG/'NS^^I 
COMMON /AFB2IN/ ALP, DELBETA ^ 

COMMON /C0NSTANT1 / SPEED, OMEGA 
COMMON /C0NSTANT2/ ALPHA, DTOT 
COMMON / CONSTANT3/RL0MDA , Q 
COMMON /C0NSTANT4/CMPI , CK , RO , R1 
COMMON /HEIGHT/ Z1 ,Z2,Z3 
COMMON /CETA/Z1EN,Z2EN,Z3IN 
COMMON /BRANCH/BRO , BRS , BRZ , BRW 
COMPLEX* 1 6 BETA, GZ, GZZ 
C0MPLEX*1 6 H2,H21 
C0MPLEX*16 EN,'Z1EN,Z2EN,Z3EN 
C0MPI«X*16-R1 ,H11 
COMPLEX* 1 6 CK, TAU, SC I, Z IMPED 
C0MPLEX*1 6 DUM1 ,DUM2,R0,R1 ,CMPI 
COMPLEX*! 6 GBAR( 32768) 

CALL INPUT (TINF, SPIREF, RESISTANCE ) 

Q= . 00002*4 • *3 • 1 41 5926*4 . 67DO* ( 1 0 .**( SPLREF/20 . DO ) ) 

SPEED=DSQRT(1 .4DO*287.DO*TINF) 

PRINT *, 'THE FOLLOWING IS AN ECHO OF THE INPUT ' 

PRINT *, 'IN THE FOLLOWING ORDER ALP DELTA ME NP DELBETA' 

PRINT *, 'SPEED OMEGA ALPHA DTOT Z1 Z2 Z3 RESISTANCE’ 

PRINT *,'TINF,SPIREF' 

PRINT *, ' ' ! SKIP A LINE 

PRINT *, 'THESE VALUES ARE ALSO WRITTEN TO FILE 44' 

PRINT *, ALP, NS, ME, N1, DELBETA, SPEED, OMEGA, ALPHA, DTOT 
PRINT *,Z1 ,Z2,Z3, RESISTANCE, TINF.SPIREF 

WRITE (44,*) 'ECHO OF INPUT ALP NS ME DELBETA SPEED OMEGA ALPHA 
& DTOT Z1 Z2 Z3 RESISTANCE TINF SPIREF' 
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VRITE(44,*) ALP, NS, ME, N1 ,DELBETA,TINP,SPIi^EP, OMEGA, ALPHA, 

& DTOT,Z1 ,Z2,Z3, RESISTANCE, Q, SPEED 
RLEMDA=OMEGA/ ( SPEED*ALFHA ) 

PI=3. HI 5926 
CMPI=(O.O,1.D0) 

C FUNCTIONS TO COMPUTE THE COMPLEX IMPEDENCE PROM CHESSELL MODEL 
PR=0MEGA/2.D0/PI 
RATIO=FR /RESISTANCE 
R=1 . +9 • 08D0*RATI O** ( - . 75DO ) 

X=>-11 .9DO*RATIO**(-.73DO) 

ZIMPED=DCMPLX(R,X) 

BRO=OMEGA/SPEED*DSQRT ( 1 .D0/(1 .DO V ALPHA*Z1+DT0T)) 
ERZ=OMEGA/SPEED*DSQRT , (- (1 .D0+ALPHA*Z2)/(1 .DO+-ALPHA*22+DTOT) ) 
BRS=CMBGA/SPEED*DSQRT ( (,1 . D0fALPHA*Z3 ) /( 1 . DCH-ALPHA*Z3+DT0T ) ) 
BRW=OMBGA/SPEED ' ‘ 

PRINT *, • * 

PRINT *, • • / f 

PRINT *, ’THE BRANCH CUTS ARE' ,BRO,BRZ,BRS,BRW 
VRITE(44,*) 'THE BRANCH CUTS ARE' ,BRO,BRZ,BRS,BRW 
DO 1,1=1 ,N1 

BETA=DCMPLX ( DFLOAT ( 1-1 ) , (-ALP))*(DELBETA) 

IF (ABS(BETA) .LE. .00000001) THEN 
BETA=DCMPLX( . OOOOOOTDO, ( -ALP ) ) *DELEETA 
END IF ' 

CALL GZALL1 (Z1 ,BETA,GZ,GZZ,JN) 

Z1HJ=M 
DUM1 =GZ 
DUM2=GZZ 4 

CALL GZALL1(Z2,BETA,GZ,GZZ,M) 

Z2M=m 

CK=Q/12.D0/CMPI/(RLEMDA**(2.D0/3.D0))*1 .0/(CDSQRT(GZ)) 

CALL GZALL1 (Z3,BETA,GZ,GZZ,EN) 

Z5EN=M 

CK=CK*1 .DO/(CDSQRT(GZ)) 

TAU=ALPHA*RLEMDA-CMPI / 2 . D0*ZIMPED*DUM2 /DUM1 

SCI 4 IMPED# ( (3-/2. )**(2.D0/3»D0))*( (RLEMDA)**(2.DO/3.DO) )*DUM1 

CALL HALL(Z1 EN,H2,H21 ,H1 ,H1 1 ) 

DUM1 =IAU*H1 +CMPI*SCI*H1 1 
DUM2=TAU*H2+CMPI*SCI*H21 
R0=-DUM1/DUM2 
DUM1 =CMPI*PI/6.D0 

R1=(CDEXP(-DUM1 )*(CMPI*RO) )/( (CDEXP(DUM1 )**2.D0)+(R0**2.D0) ) 
CALL REGION(BETA, IRBGION) 

GO TO(10,20,30,40,50,60,70,80),IREGION 
C * ** «***«** region 1 BEGINS HERE FOR Z2>Z3 OR Z>S ********** 

10 CONTINUE 

CALL HALL2(Z3EN,H2,H1 ) 

DUM1=H1+R0*H2 

CALL HALL2(Z2EN,H2,H1 ) 

GBAR ( I ) =CK*H2*DUM1 
GOTO 500 

C ««mm REGION 2 BEGINS HERE FOR Z2>Z3 OR Z>S ********** 



20 CONTINUE 

CALL HALL2(Z3M,H2,H1 ) 

DUM1=H1+R1*H2 

CALL HALL2(Z2EN,H2,H1 ) 

GEAR ( I ) =CK*H2*DUM1 
GOTO 500 

C ** * ******* REGION 3 BEGINS HERE FOR Z2>Z3 OR Z>S 
30 CONTINUE 

CALL HALL2(Z3EN,H2,H1 ) 

DUM1 = (HI +R0*H2)*R1 *CDEXP(CMPI*PI/3.D0) 

CALL HALL2(Z2M,H2,H1 ) 

GEAR ( I ) =CK*H2*DUM1 
GOTO 500 - % >. 

C «»"« » "» »»» REGION 4' BEGINS HERE FOR Z2>Z3 OR Z>S ********** 
40 CONTINUE • \ 

CALL HALL2(Z3M,H2 f H1 ) 

DUM1 =CDEXP ( CMPI*Pl/3 -'DO ) *R1 * (HI +R0*H2 ) 

CALL HALL2(Z2HJ,H2,H1 ) 

GEAR ( I ) =CK*DUMi * ( H2+ ( CDEXP ( CMPI *PI /3 • DO ) *H1 ) ) 

GOTO 500 

C ******** ** REGION 5 BEGINS HERE FOR Z2<Z3 OR Z<S ********** 
50 CONTINUE 

CALL HALL2(Z2ffl,H2,Ht) 

DUM1=Hl4B0*H2” 

CALL HALL2(Z3IN,H2,H1 ) 

GBAR ( I ) =CK*H2*DUM1 
GOTO 500 

C ******** * * REGION (, BEGINS HERE FOR Z2<Z3 OR Z<S ********* * 
60 CONTINUE 

CALL HALL2(Z2EN,H2,H1 ) 

DUM1 =H1 +R1 *H2 

CALL HALL2(Z3EN,H2,H1 ) 

GBAR ( I ) =CK*H2*DUM1 
GOTO 500 

C ********** REGION '7 BEGINS HERE FOR Z2<Z3 CR Z<S ********** 
70 CONiaiUE 

DUM1 =CDEXP(CMPI*PI/3.D0)*R1 
CALL HALL2(Z2M,H2,H1 ) 

DUM1 =DUM1 *(H1 +R0*H2 ) 

CALL HALL2(Z3IN,H2,H1 ) 

GBAR ( I ) =CK*H2*DUM1 
GOTO 500 

C ********** REGION 8 BEGINS HERE FOR Z2<Z3 OR Z<S * ** * ****** 
80 CONTINUE 

DUM1 =CDEXP(CMPI*PI/3.D0)*R1 
CALL HALL2(Z2EN,H2,H1 ) 

DUM1=DUM1*(H1+R0*H2) 

CALL HALL2(Z3EN,H2,H1 ) 
GBAR(I)=CK*(H1*CDEXP(CMPI*PI/3.D0)+H2)*DUM1 
GOTO 500 

C HTO OF GBAR (BETA) CALCULATIONS BASED ON REGIONS DETERMINED 

C MULTIPLY BY BETA ONLY TO MATCH VAN MOORHEM DEFINITION TO 
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C RICHARDS AND ATTENBOROUGH. 

500 CONTINUE 

GBAR(I)=GBAR(I)*BETA 

1 CONTINUE 

PRINT V LAST BETA EXCUTED IS '.BETA 
CALL DAFB2(GBAR) 

DO 2,1=1 ,Nl/2 

RAD=2 . DO*PI *DPLOAT (1-1 )/(DPL0AT(N1 )*DELHETA) 
IP(RAD.LE.O.O ) THEN 
GOTO 5 
END IP 

RAD2=DL0G1 O(RAD) 

G=k>.DO*DL0G10(ABS(£BAR(I))) V 
WRITE (44, 9095) RAD2,G 
5 CONTINUE 

2 CONTINUE 

PRINT *, 'FINALLY FINISHED' 

9093 F0RMAT(5X,3G18.8) f 

STOP », 

END 


SUBROUTINE INPUT (TINF,SPLREF,RESISTANCE) 

Q »»** » ■ « KKXXKKKKKKKKKKKKXKKKKKKKXK X * « «-* « «-* « * « iHHHHHHHHH H Ht HHH H HHHHHK < H HHt 

C THE PURPOSE OF THIS ROUTINE IS TO INPUT ALL NECESSARY PARAMETERS TO * 
C THE MAIN ROUTINE. THE DEFINITION OF EACH PARAMETER WILL BE DEFINED * 
C AT ITS RESPECTIVE INPUT*. * 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

common /mem/ ns.me.ni 

COMMON /AFB2IN/ ALP.DELBETA 
COMMON /C0NSTANT1 / SPEED, OMEGA 
COMMON /C0NSTANT2/ ALPHA, DTOT 
COMMON /C0NSTANT3/ RLEMDA.Q 
COMMON /C0NSTANT4/. CMPI,CK,R0,R1 
COMMON /HEIGHT/ El ,Z2,Z3 
COMMON /CETA/ Z1M,Z2EN,Z3EN 
COMMON /BRANCH/ BRO,ERS,BRZ,BRW 
WRITE (6,899) 

WRITE(6,900) 

READ *, TINF 
WRITE (6,901) 

READ *, DTOT 
WRITE (6,902) 

READ *, ALPHA 
WRITE (6,903) 

READ *, SPLREF 
WRITE (6,904) 

READ *, OMEGA 
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900 


901 

902 

903 

904 

905 

906 

907 


906 

909 

910 
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WRITE (6,905) 

READ *, RESISTANCE 
WRITE (6,906) 

READ *, Z2,Z3 
Z1=0.D0 
WRITE (6,899) 

WRITE (6,907) 

READ *, ALP 

\ ! PESSARY DUE TO PREVIOUS DEFINITION USED 
WRITE \Oj906J 
READ *, NS 
WRITE (6,909) 

READ *, ME * v 

WRITE (6,910) . 

READ *, N1 

WRITE (6,911) ‘ % 

READ *, DELBETA . .. 

•///////////////) .! CLEARS THE SCREE® 

P0RMT(’0',///,l6x,3eH THE FOLLWINO ARE PHYSICAL VARIABLES. 

,//,12X,43H INPUT THE TEMPERATURE AT Z EQUAL INFINITY. 

ft) 

FORMAT (XJ2X,45H INPUT THE TEMPERATURE CHANGE DELTA T/T(INF) . 

FORMAT (X^12X,46H INPUT THE TEMPERATURE PROFILE CONSTANT ALPHA. 

FORMAT (X^12X,46H INPUT THE REFERENCE SOUND PRESSURE LEVEL. 

FORMAT (XJ2X.46H INPUT THE ANGULAR FREQUENCY IN RADIANS /SEC . 

FORMAT (XJ2X,47H INPUT THE CHESSELL MODEL FLOW RESISTANCE (cgs) 

FORMAT (X,12X,46H INPUT RECEIVER AND SOURCE HEIGHTS SEPARATED 
,/,25H BY A COMMA.,/) 

FORMAT (///,10X,37H THE FOLLOWING ARE NUMERIC VARIABLES. 

,//,12X,39H INPUT ALP. THE IMAGINARY PART OF BETA. 

• <1 ,/,12X,43H NOTE POSITIVE VALUES ARE ABOVE THE AXIS. /) 
FORMAT (X,T2X,46H INPUT THE DELTA PARAMETER USED IN INTEGRATION 

ft) 

FORMAT (X,12X,46H INPUT THE NUMBER OF TERMS IN SUM TO INFINITY. 

FORMAT (XJ2X,46H INPUT THE NUMBER OF POINTS TO BE USED. 

FORMAT (X,12X,46H INPUT THE STEP SIZE IN BETA (OR K) SPACE. 

RETURN 

END 


SUBROUTINE REGION (BETA , IREGION ) 


c 

c 


SUBROUTINE TO DETERMINE WHICH OP 8 REGIONS 
THE EVALUATION IS TO TAKE PUCE IN. 


* 

* 


C 

c 

c 

c 

c 


c 


IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

INTEGER DOUBLE PRECISION K1 
C0MPLEX*16 BETA.G32C1 ,G32C2,G32C3 
C0MFLEX*16 G32PI,G32PI2,G32PI3 
COMMON /INTEG/NS,ME,N1 
COMMON /AFB2IN/ALP,DELBETA 
COMMON /CONSTANTI /SPEED, OMEGA 
COMMON / C0NSTANT2/ ALPHA , DTOT . 

COMMON / CONSTANT3/RLEMDA , Q v 

COMMON /C0NSTANT4/CMP1 ,CK,RO,R1 
COMMON /HEIGHT/ Z1 , Z2 , Z3 
COMMON /CETA/ZtEN,Z2EN,Z3EN 
COMMON /BRANCH/ERO* BRS r BRZ , ERV * 

INTEGER IREGION 
CALL G32ALL(Z1 ,fcETA*G32C1 ) 

CALL G32ALL ( Z2 , BETA , G32C2 ) 

CALL G32ALL(Z3,BEIA,G32C3) 

IN THE EVENT THAT BRZ OR BRS AND ERW ARE VERY CLOSE THE FOLLOWING 
IS USED TO INSURE PROPER , REGION IS CHOSEN. VI AND V2 ARE SIMPLE 
WEI SIT' FACTORS TO, DETERMINE WHICH DIRECTION TO EVALUATE REGIONS. 

BY WEIGHTING THE SELECTION DEPENDING HOW CLOSE BETA IS TO BRW 
THE REGIONS ARE FOUND IN THE REVERSE ORDER. 

POSMAX=MAX(Z2,Z3) 

¥ 1 * 1.0 , 

V2=1 .0 * 

DUD=OMEGA/SPEEB • 

DUD=DUD*W1 +DUD*W2*DSQRT ( ( 1 . 0+ALPHA*P0SMAX) / ( 1 +ALPHA*POSMAX+DTOT ) ) 
DUD=DUD/2 
IF(Z2.LT.Z3) THEN 
GOTO 120 
END IF 

FOR **(Z2>Z3 OR Z>S)** THE FOLLOWING IDENTIFY REGIONS OF SPACE 
All =DIMAG(G32C1 ) 

AI2=DIMAG(G32C3) 

AI3=DIMAG(G32C2) 

IF(REAL(BETA) .gt.dud) then 
GOTO 110 
END IF 

IF(AI3-GT.0.0) THEN 
IREGI0N=4 
GOTO 150 

IF(AI2.GT.0.0) THEN 
IREGI0N=3 
GOTO 150 
END IP 

IF(AII.GT.O.O) THEN 
IREGI0N=2 
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end ip 

IREGI0N=1 
GOTO 150 
110 CONTINUE 

IF(AI1 .IT. 0.0) THEN 
IREGIQN=1 
GOTO .150 
mu IP 

IP(AI2.LT.O.O) THEN 
IREGIQN=2 
GOTO 150 

END IP • <; 

IPCAIJ.Iff.O.O) THEN ’• 

IREGI0N=3 
GOTO 150 

END IP . , . 

IREGI0N=4 

GOTO 1 50 ■ } , 

C F0R**(Z2<Z3 OR Z<S)** THE FOLLOWING IDENTIFY REGIONS OP SPACE 
1 20 CONTINUE 

All =DIMAG(G32C1 ) 

AI2=DIMAG(G32C2) 

AI3=DIMAG(G32C3) . 

IP (REAL( BETA ).GT. DUD) THEN 
GOTO 130 
END IP 

IP (AI3.GT.0.0) THEN 
IKEGI0N=8 
GOTO 150 
END IP 

IP (AI2.GT.0.0) THEN 
IREGI0N=7 
GOTO 150 
END IP 

IP^ffll .GTtO.O) THEN 
IREGI0N=6 
GOTO 150 
END IP 
IREGI0N=5 
GOTO 150 
130 CONTINUE 

IP (All .LT.O.O) THEN 
IREGI0N=5 
GOTO 150 
END IP 

IP (AI2.LT. 0.0) THEN 
IREGI0N=6 
GOTO 150 


IREGI0N=7 
GOTO 150 
END IP 
IREGI0N=8 
GOTO 150 
150 CONTINUE 
RETURN 
END 


C 

C 

C 

C 


SUBROUTINE G32ALL(Z,BETA,G32C) 

g KKKKKKK IH H( > > K K K K ) HHK -> XllKKK ^ KKKKKKKItKKKK 


V 

OHHnmnmnnKH 


c* • * 

C* G32ALL ‘CALCULATES g>/2. FUNCTION * 

C* ■ t * 

C* \ * 

C<H HKK < Um H HHHmi(KK * «« l lKltlCKK * «K«>(KKKIH t KKKKKKKKIHUKKK)lKKK -) HHm(«)m 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

C0MPLEX*16 BETA, PHI, SORT 1 ,SQRT2,F00,G32C,M0DL0G 

C0MPLEX*1 6 SI ,A,B,S 

COMMON /CONSTANT 1 /SPEED, OMEGA 

COMMON / C0NSTANT2 / ALPHA , DTOT 

COMMON /BRANCH/BRO , BRS , BRZ , BRW 

AA=1 .+ALPHA«2+DT0T 

B=1 .-( ( SPEED*EETA/OMEGA ) * ( SPEED*BETA / OMEGA ) ) 

SI =SQRT1 (BETA,Z) 

S=SQRT2(EETA) 

PHI=S1 /(S*DSQRT(AA)' ) 

IF(DTOT.EQ..O.OR.CDABS(1 .-PHI). LT. ID-8) THEN 
F00=.0 


ELSE 

P00=M0DL0G((1.+PHI)/(1 .-PHI)) 

ENDIF 

G32C=#SQRT CAA )*S1 -. 5*DT0T*F00/S 

RETURN 

END 


C 

C 

c 

C 

c 

c 

c 

c 


BEGIN OF FUNCTIONS USED ABOVE. 


SQRT1 AND SQRT2 ARE FUNCTIONS TO CALCULATE SORT (BETA) GIVEN 
DESIRED BRANCH CUTS AND DIRECTION +IMAGINARY OR -IMAGINARY. 


* 

* 


FUNCTION SQRT1 (BETA,Z) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

C0MPLEX*1 6 BETA.AA.SQRTI 

COMMON /C0NSTANT1 /SPEED, OMEGA 

COMMON /CONST ANT 2 /ALPHA , DTOT 

AA=1 .-( ( (SPEED/OMBGA)*BETA)*( ( SPEED/OMEGA )*BETA) ) 
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BB=1 . +ALPHA*Z+DTOT 

BRAN1 = (OMBGA/SPEED )*DSQRT ( ( 1 . +ALPHA*Z ) / ( 1 . +ALPHA*Z+DTOT ) ) 
IF ((DREAL(BETA) .GE. BRAN1 ) .AND. 

1 (DIMAG(BETA) .IE. 0.0)) THEN 
SQRT1 =-CDSQRT ( AA*BB-DT0T ) 

ELSE 

SQRT1 =CDSQRT(AA*BB-DTOT) 

END IP 
RETURN 
BID 




FUNCTION SQRT2(EETA) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

C0MPLEX*1 6 BETA,AX,SQRT2 
COMMON /BRANCH/BRO , BRS , ERZ , BRW 
COMMON /C0NSTANT1 /SPEED, OMEGA 
COMMON /C0NSTANT2 /ALPHA , DTOT 
AA=1 .-( ( ( SPEED/OMEGA )*BETA)*( ( SPEED/OMEGA )*BETA) ) 
IP ( (DREAL(BETA) .GE. BRW) .AND. 

1 (DIMAG (BETA) .12!.. 0).) THEN 
SQRT2=-CDSQRT ( AA ) 

ELSE 

SQRT2=CDSQRT ( AA ) 

END IP 
RETURN 

END ' 


FUNCTION MQPLOG COMPUTES THE LOG OP BETA GIVEN DIRECTION AND * 

LOCATION OP BRANCH CUTS. * 

FUNCTION MODLOG(QUAN) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

C0MPLEX*16 QUAN,MODLOG 

IP ((DREAL(QUAN) .LE. 0.0) .AND. (DIMAG(QUAN) 

1 .GE. 0.0)) THEN 

M0DL0G=CDL0G ( QUAN ) +DCMPLX (0 . 0 , - 2 * 3 .141 5927 ) 

ELSE 

MODLOG=CDLOG ( QUAN ) 

END IP 
RETURN 
END 


END OP FUNCTIONS USED ABOVE. 


SUBROUTINE HALL(Z,H2,H21 ,H1 ,H11 ) 

HALL USES SUBROUTINE CGBAIR TO CALCULATE 1/3 ORDER 
HANKEL FUNCTIONS FROM AIRY FUNCTIONS. 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMPLEX*! 6 Z,AI,BI,AIP,BIP,K,KS,H1 ,H2,H11 ,H21 
C0MPLEX*16 ARG,CI 
CI= DCMPLX(O.DO,1.DO) 

PI= 3*141 5926 54D0 

ARG= DCMPLX(O.DO,-PI/6.DO) \ 

K= (12.D0)**(1.D0/6.D0 , )*CDEXP(ARG) 

KS= DCONJG(K) * ' • 

CALL CGBAIR (-Z,AI,BI,AIP;BIP) 

H1= K*(AI-CI*BlJ 
H2= KS*(AHCI*BI), : * 

H1 1 = -K*(AIP-CI*BIP) 

H21 = -KS*(AIP4CI*BIP) V 

RETURN 

END 


SUBROUTINE CGBAIR (Z,AI,BI,AIP7BIP) 

<hhh h hhhhhhhhhhh h hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhi 

CALCULATE AIRY FUNCTIONS FOR COMPLEX*! 6 ARGUMENT * 

REF. HANDBOOK OF MATHEMATICAL FUNCTIONS, AERAMOWITZ AND STEGUN. * 
ENTRY: : * 

CALCULATE ARGUMENT(Z) AND ABSOLUTE VALUE(Z) * 

IF /Z/ LT 6 * 

THEN USE EQS. 10.4-2 THRU 10-4-5 FOR AI,BI,AIP,BIP * 

10 EISE IF ARG(Z) LT Pl/3 * 

THEN CALCULATE ZETA(Z) * 

USE EQS« 10.4.59, 10.4.61, 10.4.63, 10.4*66 FOR AI,BI,AIP,BIP* 

20 ELSE CALCULATE ZETA(-Z) * 

USE EQS. 10.4.60, 10.4.62, 10.4*64, 10.4.67 FOR AI,BI,AIP,BIP* 

ENDIF * 

ENDIF * 

EXIT * 

END * 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

C0MPLEX*16 Z,AI,BI,AIP,BIP,ZETA,CZETA,Z14,SUM1 ,SUM2,SUM3,SUM4, 
1 ZETAP,FACT1 ,FACT2,SN,CS,FTERM,FPTERM,GTERM,GPTERM,F,FP,G,GP,Z3 
C0MPLEX*16 VZETA.VZETAP 
DIMENSION C(21),D(21) 

DATA Cl ,C2,PIRT,PI4/.355O280539DO, .25881 94038D0,1 .772453851D0, 
+ *7853981 635DO/ 

+ *0371 33487654321 DO, *0379930591 27800D0, 
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1 . 0576491 9041 2669DO, . 1 1 609906402551 DO, 

+ .291 591 39923O74D0, .87766696950998DO, 

2 3*0794550301 731 DO, 1 2. 341 573332345DO, 

+ 55 • 62278536591 4D0, 278. 46508077759DO, 

3 1533*1 6943201 27D0, 9207 * 2065997258D0 , 

+ 59892. 51 3565875D0, 41 9524*8751 1 653DO, 

4 31 48257* 41 78666D0, 251 9891 9*871 601 DO, 

+ 21 4288036. 96366D0, 1 929375549- 1 823DO, 

5 1 8335766937 • 889D0/ 

DATA D/1 .DO, 

+ -.097222222222221 DO,-. 0438850308641 97DO,-.O42462830789894D0, 

1 -.0626621 63492031 DO, -. 1 241 0589602727D0, -. 308253764901 OTDO, 

2 92047999241 291 DO, -3 .€1 049358464851)0,-1 2 .807293O6O735DO, 

3 -57 . 50830(351 391 1 DO, -287 . 0332371 0920D0, -1 576 . 3573033370DO, 

4 -9446 . 354823O953D0 , -6t 335 * 706663847D0, -428952 . 40040004D0 , 

5 -321 4536. 521 4006D0, -25697908. 383909DO, -21 8293420.8321 4D0, 

6 -1 963523788. 9909D0, r 1 8643931 088 . 1 05D0/ 

ABSZ=ABS(Z) 

IF(ABSZ.BQ.O) GO $0 5- 

IP(ABS(DIMAG(Z)).LE.1 .D-12.AND.DREAL(Z).IZT.0.D0) GO TO 5 
ARGZ=ATAN2 (DIMAG(Z ) , DREAL( Z ) ) 

GO TO 4 

3 ARGZ=O.DO 
GO TO 4 

5 ARGZ=3 . 1 41 5926535898D0 

4 CONTINUE 

IF(ABSZ.GT.6.D0) GO TO 10 
ASCENDING SIRIEg 
EQS. 10.4.2,10.4.3 
CONTINUE '* 

Z3=Z**3 

FTERM=1 .DO 

FPTERM=Z*Z/ 2 . DO 

GTEEM=Z 

GPTERM=1 .DO 

GMM=1.D-13*ABSZ 

F=FTERM 

FP=FPTIRM 

G=GTERM 

GP=GPTERM 

KKKT=100 • ! ADJUST KKKT TO INSURE CONVERGENCE IP NECESSARY 
DO 1 1=1, KKKT 
13=3*1 

FTERM=FTERM*Z3/((I3-1 .D0)*I3) 

PPTIRM=PPTERM*Z3/( 13* ( 13+2 . DO ) ) 

GTERM=GTERM*Z3/( I3*( 13+1 .DO)) 

GPTERM=GPTERM*Z3/ ( ( 13-2 . DO ) *13 ) 

F=F+FTERM 

PP=PP+PPTERM 

G=G+GTERM 

GP=GP+GPTIRM 

IP(ABS(GTERM) .LE.GLIM) GO TO 2 


O O Cl o ooo 


1 CONTINUE 
PRINT 6000, Z 

6000 F0RMAT(/’ Z='2E14.5, ' ERROR IN CGBAIR, NONCONVERGENCE ' ) 

2 AI=C1*F-C2*G 
AIP=C1 *FP-C2*GP 

BI=1 .732O5O8O8D0*(C1*P4C2*G) 

BIP=1 .752050608D0*(C1 *FP+C2*GP) 

GO TO 9999 

C ASYMPTOTIC EXPANSIONS FOR /Z / LARGE 

10 SIGN=1 .DO 
SUM1=0.D0 
SUM2=0.D0 

SUM5=0.D0 * v 

SUM4=0.D0 

PIBY3=3 • 1 41 5926D0/3 .DO-- 
IF(ABS(ARGZ).GE:PIBY3) GO TO 20 
/ARG(Z)/ LE PI/3 • t 
EOS. 10.4.59, 10.4.61, 10.4.63, 10.4.66 

ZETA=CZETA ( ABSZ , ARGZ ) * 

DO 11 1=1,12 
K=I-1 

ZETAP=ZETA**K 
SUM1 kSUMI +SIGN*C(I ) /ZETAP 
SUM2=SUM2+SIGN*D( I ) /ZETAP 
SUM3=SUM3+C(I)/ZETAP 
SUM4=SUM4+D(I ) /ZETAP 

11 SIGN=-SIGN 

Z1 4=ABSZ** . 25DC^DCMPLX ( COS (ARGZ/4 .DO) , SIN(ARGZ/4.D0) ) 
FACT1 = . 5D0*EXP ( -ZETA ) / ( PIRT*Z1 4 ) 

FACT2=. 5D0*EXP(-ZETA)*Z1 4/PIRT 
AI=FACT1 *SUM1 
AIP=-FACT2*SUM2 
FACT1 =EXP ( ZETA ) / ( PIRT*Z 1 4 ) 

FACT2=EXP ( ZETA )*Z1 4/PIRT 
BI=FACT1 *SUM3 
BIP=FACT2*S0M4 
GO TO 9999 

/ARG(Z)/ GT PI/3 NOTE CHANGE ABOVE 
EQS. 10.4.-60, 10.4.62, 10.4.64, 10.4.67 

20 CONTINUE 

ARGZ=ATAN2 (-DIMAG(Z ) , -DREAL(Z) ) 

ZETA=CZETA ( ABSZ , ARGZ ) 

VZEEA=1 .DO/ZETA 
LLL=10 

DO 21 1=1 ,LLL 
K 2 = ( 1 - 1 )*2 
J=K2+1 

VZETAP=VZETA**K2 
SUM1 =SUM1 +SIGN*C( J )*VZETAP 
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SUM2=SUM2+(SIGN*C(J+1 )*VZETAP*VZETA) 
SUM3=SUM3f(SIGW*D( J )*VZETAP) 

SUM4=SUM4+ ( SIGN*D ( J+1 )*VZETAP*VZETA) 

21 SIGW=-SIGN 

Z1 4=ABSZ** . 25D0*DCMPLX ( COS ( ARGZ/ 4 • DO) , SIN ( ARGZ /4 • DO) ) 
FACT1 =1 . DO/ ( PIRT*Z1 4 ) 

FACT2=Z1 4/PIRT 
SN=SIN ( ZETA+PI4 ) 

CS=C0S (ZETA+PI4) 

AI=FACT1 *(SN*SUM1 -CS*SUM2) 

AIP=-FACT2*( CS'*SUM3+SN*SUM4) 

BI=FA'CT1 *(CS*SUM1 +SN*SUM2) X 

BIP=FACT2*( SN*SU{6-CS*kJM4 ) 

9999 RETURN 

END . „ X 

v * 

; > • 

BEGIN OF FUNCTIONS USED ABOVE 


FUNCTION CZETA (ABSZ , ARGZ ) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMPLEX* 1 6 CZETA 
ARG=ARGZ*1 .5DO * 

CZETA= ( ABSZ** 1 .5DO)*DCMELX(COS(ARG) ,SIN(ARG) )*.66666666666667DO 

RETURN 

IND 


MD OF FUNCTIONS USE) ABOVE. 


SUBROUTINE GZALL1 (Z,BETA,GZ,GZZ,IN) 

Q X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X 

C GZALL1 CALCULATES ALL THE PARTIAL DERIVATIVES * 

C OF THE g FUNCTIO N. * 

IMPLICIT DOUBLE PRECISION(A-H,0-Z) 

COMPLEX *16 BETA,GZ,GZZ,G,GB,GEB,SQRT1,IN 
COMMON /CONSTANT 1 /SPEED, OMEGA 
COMMON /CONSTANT2/ALPHA,DTOT 
CALL GALL(Z,BETA,G) 

A=1 . +ALPHA*Z+DTOT 

BRAN =OMEGA /SPEED*SQRT ( ( A-DTOT ) /A ) 

IF (DREAL(G).LE. .O.AND.DIMAG(G) .GE.O. ) THEN 
SI=-1 • 

ELSE 

SI=1. 

hlNl) J j* 1 

GZ=SI*2 . *ALPHA*SQRT 1 (BETA, Z ) / (3 • *CDSQRT (G*A) ) 

C=2 . * ALPHA **3 • D0*DT0T/ ( 9*A**2 . DO ) 

G ZZ=C/(GZ*G)— 5*GZ**2. DO/G 


oooo ooooo oooo o o o oooo 


T= ( 3 . *0MBGA/ ( 2 . *ALPHA*SPEED ) ) ** ( 2 . DO/3 . DO ) 

EN=T*G 

RETURN 

END 


SUBROUTINE GALL (Z, BETA, G) 

GALL EVALUATES THE g FUNCTION * 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMPLEX *16 G,BETA,G32 
COMMON /CONSTANT 1 /SPEEL,- OMEGA 
COMMON / CONST AOT2 /ALPHA , DTOT 
CALL G32ALL(Z,BETA # G32)- 
G=CDEXP( 2 . /3« *CDI)0G(G32 ) ) 

RETURN • | .. 

END 


SUBROUTINE HALL2(Z,H2,H1 ) 

KKKXKKKK K XKKKKKKKKKK iHHHH H HHHHHt X-K K MHK KHHHt IHHHHHHHHHHHHHHHHHt 

HALL2 USES SUBROUTINE CGBAIR TO CALCULATE 1/3 ORDER * 

HANKEL FUNCTIONS EROM AIRY FUNCTIONS. NOT THE * 

DERIVATIVES AS HALE DOES. * 

IMPLICIT DOUBLE ’PRECISION (A-H.O-Z) 

C0MPLEX*16 Z, AI,BI,AIP,BIP,K,KS,H1 ,H2,H1 1 ,H21 
COMPLEX*! 6 ARG,CI 
C0MPLEX*16 BETA 
CI= DCMPLX(O.DO,1.DO) 

PI= 3.^5926541)0 

ARG= -‘DCMPLXCO: DO, -Pl/6. DO ) 

K= (12.D0)**(1.D0/6.D0)*CDEXP(ARG) 

KS= DCONJG(K) 

CALL CGBAIR (-Z,AI,BI,AIP,BIP) 

HI = K*(AI-CI*BI) 

H2= KS*(AI+CI*BI) 

RETURN 

END 


SUBROUTINE DAFB2(F) 

C UWyM VWyyMMy|f>fWWMMM>(y|f WM«WHM«WVI/MWVM»V I( «VVIfyMV1tW^MI>M>fMWMM«WXWWW>/WWWMMM 

C SUBROUTINE TO ACCURATELY DO THE HANKEL TRANSFORM OF THE SOUND * 
C PRESSURE LEVEL. * 


C F(NP)=GBAR(NP) MUST BE SAMPLED AT NP POINTS WITH K=(N-1 ,ALP) * 
C ALP REPRESENTS THE DISTANCE ABOVE THE REAL AXIS THE FUNCTION WILL * 
C BE INTEGRATED. * 
C NS IS A PARAMETER REPRESENTING ADDITION OF AN ANALYTICAL FUNCTION * 
C TO F(NP) * 
C M IS THE NUMBER OF TERMS USED TO APPROXIMATE F(NP) TO INFINITY * 


IMPLICIT DOUBLE PRECISION (A-fl.O-Z) 

COMMON /INTBG/NS,ME,N1 
C0MM0N/AFB2IN/ALP, DELBETA 
COMPLEX*! 6 F(N1),CF,CARG,SUM,FNP,CMPI,D1 
NP=N1 \ 

CMPI=DCMPLX(O.DO,1 .DO) 

DELK=DELEETA ‘ • 

PI=5*141 5926DO „ '' 

C SUBTRACT THE ANALYTICAL FUNCTION IF- NS > ZERO 
C ADJUSTING THE SUBTRACTION MULTIPLIER CF 
IF(NS.LE.O) GOTO >11 
CF=DCMPLX(0.D0,0.iX)) * 

IF (ALP.EQ.O.O) THEN 
CF=DFLOAT (NP) /DFLOAT (NS ) 

CF=CF*F(2) 

END- IF 

IF (ALP.NE.O.O) THEN 

CF=CMPI*DFLOAT(NP)*F( 1 )/(DFK3AT(NS)*ALP) 

END IF 

C SUBTRACT THE ANALYTICAL FUNCTION IF NS>0 
DO 10,1=1 ,NP » 

D1 =DCMPLX ( DFLOAT ( 1-1 ) , (-ALP) ) 

CARG=DFL0AT ( NS) * ( -D1 ) /DFLOAT (NP) 
F(I)=F(I)-CF*(1 .DO-CDEXP(CARG)) 

10 CONTINUE 

1 1 CONTINUE 

IF(ALP.EQ.O.O) F(1 )=DCMPLX(0.D0,0.D0) 

FNP=F(NP) 

DO 12fI=2,NP. 

D1 =DCMPLX ( DFLOAT (1-1 ) , (-ALP) ) 
F(I)=F(I)/(CDSQRT(D1 )) 

12 CONTINUE 

IF(ALP.NE.O.O) F(1 )=F(1 )/CDSQRT((-CMPI)*(ALP)) 
C ADD TERMS TO INFINITY IF ME>0 
IF(ME.IT.1 ) GOTO 20 
DO 15,1=1 ,NP 

D1 =DCMPLX ( DFLOAT (1-1 ) , (-ALP) ) 

CF=D1 /DFLOAT (NP) 

CALL ZETA(NP,ME,CF,SUM) 

F(I)=F(I)+FNP*SUM 
15 CONTINUE 
20 CONTINUE 
C DO THE FFT 

CALL F0RK(NP,F,1) 

C ADD ALTERNATE TERMS TO GIVE NP/2 SAMPLES TRANSFORMED 
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CF=DEUC*DPLOAT (NP )* ( CDSQRT ( -CMPI ) )/(2.D0*PI ) 
DO 25,I=2,NP/2 

A=DEXP ( DPLOAT ( 1-1 ) * ( ALP ) * 2 . D0*PI /DPLOAT (NP ) ) 
P(I )=A*P(I )4CMPI*P(NP-I+2)/A 
F(I)=F(I)*CF/DSQRT(DFL0AT(I-1 )) 

CONTINUE 

RETURN 

END 


C 
C 
C 

C ' ' 

SUBROUTINE FORK(IX,CX,’SIGNI) . 

C «««»«« mUKO « Kll > KKKmCKK«IUK I H ( JH« » m»,yKJm«l'KKK«KKKKKKKKmKKKK MICK mm m m 

C A PAST FFT GIVEN BY J.F.. CIAERBOUT, "FUNDAMENTALS OP GEOPHYSICAL * 
C DATA PROCESSING" . * 

Cunmimnimmi »» muT < » < x x x k k x m> x kx x k x k « > « m k x k kk » x x » x x k k k x k x »k 
IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

C0MPLEX*16 CX(LX),CARG,CW,CTEMP,CI,DUM1,CMPI 

INTEGER SIOTt 

J=1 

CMPI=DCMPLX(O.DO, 1 .DO) 

PI=3. HI 5926 . . ' 

SC=DSQRT(1 . DO/DPLOAT ( LX ) ) 

DO 30,1=1, LX 
IP(I.GT.J) GOTO 10 
CTEMP=CX(J)*SC 

cx(j)=cx(i)*sd 

CX(I)=CTEMP - 
10 M=LX/2 

20 IP(J.LE.M) GOTO 30 

J=J-M 
M=M/2 

IP(M.GE.I) GOTO 20 
30 . . J=J4M 

L=1 4 ’ 

40 ISTEP=2*L 

DO 50,M=1 ,L 

CARG=CMPI*PI*DBLE(SIGNI )*DBLE( (M-1 ) ) /DBLE(L) 

CW=CDEXP (CARG ) 

DO 50,I=M,LX,ISTEP 
CTEMP=CW*0X(I+L) 

CX(I+L)=CX(I )-CTEMP 
50 CX(I)=CX(I)+CTEMP 

L=ISTEP 

IP(L.Lr.LX) GOTO 40 

RETURN 

END 


C 

C 

C 

C 
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SUBROUTINE ZETA(NP,M,A,SUM) 

SUBROUTINE TO ADD THE NECESSARY TERMS TO 
EXPRESS INVUtTIBLE FUNCTION TO INFINITY. 

WILL USE DOUBLE PRECISION. 

SUM=SUM OF 1/(NP\5)*1/((J+A)\5) FOR J=1 TO 
INFINITY MINUS SOME CONSTANT WHICH IS 
INDEPENDENT OF A. 

tHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH H HHt 

IMPLICIT DOUBLE PRECISION (A-H f O-Z) 

COMPLEX*! 6 A,SUM,D2 
D2=DCMPLX(DFL0AT (M) , 0. DO ) , 

SUM=1.D0/(M+A) * ' ' 

SUM=2.DO*(DSQRT(DFH5AT{M))-1 .DO/CDSQRT(SUM)) 

$ -0.5*CDSQRT(SUM)*(1 .04SUM*(1 .0/1 2.0-SUM*SUM/1 92.0)) 

DO 10,J=1,M 

SUM=SUM+1 .DO/CDSQBT(J+A) 

10 CONTINUE 

sum=sum/dsqrt ( dfDoat (N? ) ) 

RETURN 

END 


APPENDIX B - The Inversion Case Program 


The code for the lapse case follows. This code is extremely similar to the lapse 
case and the subroutine names and functions, variable names and hints are identical 
or at least very similar to those in the lapse case. 


non non non n n oo nnnnr>nnnnnonnnnnnn 


PROGRAM MAIN 


ORIGINAL PAGE v ? 
OF POOR Q 7 JA uTY 


Main prograa: will call the subroutines tot 

1) Input the prograa pirantttri, 

2) Fine the points in tha g*3/2 function ahere 
region changes occur, 

3) 3uild the matrix of values obtained by aarchlng 
slightly above the real axis, 

4) Parlor* the Fankel Transfora on tha aatrix 
obteined in step three. 

5) Print the results. 


********************** 


IMPLICIT CCU5LE FRECISION CA-H/O-Z) 

COMMON /«a:n/ F 

CQMMCh /CONST/ Cl, PI 

CCMMCN /C ON ST 1 / S P E EC , OMEGA 

COMMON / CON ST2 / ALPHA/OTCT 

COMMON / C ON ST 3 / Z,ZREF,S 

COMMON / REGION/ R ST /R S2/RZT ,RZ2,R01 ,R02 

COMMON / S T A T : / FCE, RESISTANCE 

CCMMCN / I N T EG / NS/ME/NQPTS 

COMMON /AFS2XN/ hGHT,C€L8ETA 

CCMPLEX*1 b F(327€£)*CX/GL r ESS/ ROOT x G 3 2 

Cl = CCMPIX ( j • C G , 1 .00) 

PI = OAT AN2 (O.CO#-1 .DC) 

PRINT *,'?!= ',PI 


CALL INPUT (0£LocTAxHGhT/RT1/RT2) 
ANG = -PI/2.00 


PRINT ^'DETERMINING REGION CHANGE COORDINATES' 
CALL REGION. FIND UNG/RT1 ,RT2> 


PRINT *,'BUIL0y<G MATRIX' 

CALL 5UIL0MATRIX < 0 E L 3 E T A , HGHT ) 


PRINT */'DClNG * HANKEL TRANSFORM ' 
CALL HANKEL 



print *,'Prirting results to F0RC45.OAT* 
CALL PRINTALL <D6tBETA,N0PrS> 


PRINT COKPCETe**' 

STOP 
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ORIGINAL PAGE IS 

OF POOR QUALITY 



SUBROUTINE INPUT (0ELBETA/HCHT,RT1/RT2) 

****************** ** ************************************** 

I THIS IS THE INPUT SECTION OF THE PR06RAH * 

. ALL PARAMETERS NECESSARY FOR THE OPERATION * 

* Or The PROGRAM HILL BE INPUT IN THIS ROUTINE • 

* ni ,»******••*< . ... ,»*i,*********t ********** ***************** 


SOC 


implicit double 

COMMON /CONST/ 
COMMON /CQNST1/ 
COMMCN /CONST 2/ 
COMMON /CONST3/ 
COMMON /STATE/ 
COMMCN / INT cG / 
COMPLEX«1d Cl 


PRECISION (A-H/O-Z) 
CIrPi 

SPeECrOMEGA 

JLPHA/CTCT 

Z/tREF/S 

SCE/RcSISTANCE 

KS/ME/NOPTS 


.RITE 

PRINT 

PRINT 

READ 


(6# dGC) 

*/'lN 3 LT st YSICAL PARAMETERS: 
* / # 


Catactor Haight* Sourca haight* - 
end Angular Frtquincy' 

Pk^Nl */ 'INPlT% 7ATE PARAMETERS;! a mo gradiant *Ta»paratura proTil# 

# - constant and T aap - at infinity 

print * / 

kcAG *,3T0T* ALPHA/T1NF . 

PRINT ** 'Input t h a chassall irodel P LCW RESISTANCE 

k=;d resistance 

PRINT ** 'Input Xt -9 NUMBER OF DATA PCINTS 
READ **NCPT$ 

T ST =CLOG(O p LOAT(NDPTS) ) /CLOG (2. C 01 
IrCQABSCTST-CNlNTCTST)) -GT. 1-0-8) THEN 
N3PTS-DNINT(2-CG**CNINT(TST)) 

PRINT **'InvaT*d input (i«ust fca a poaar ° f ** TtY AGAIN 
P"fNT (ax. '.NOPTS. 

GOTO 1 

PRINT */' INPUT HANKEL TRANSFORM PARAMETERS: 

Haight* 

R E AO **N$/MS*HGHT 
ZREF * 0.C0 

SPEED * 05CRT < 401 • ECO * TINF) 

RT1 - OMEGA/ SPEED 

RT2 * RT1*DSCRT<1. 00/C1. CC^OTOT)> 

C EL 3 ETA = RT2 * 1.0100 /NOPTS 

U2XTE (6*800) 

FORMAT ( # Q # *2CX*1H */// /////// lltftlllilf tltft 

PRINT *#*R«faranca hght ***IRE* - 

PRINT a/^Oatactor haighf*^7X 

***Sourca Haight * # *$ ■ -j .. ^ , 

*/ * Fr tcuancy Crad7*^/ONEGA 
★**«paad of found *%S*fEO 
*/ f T««p Gradiont * # #OTOT 


NS * ME and Intagration 


f eltirt tha acraan 


PRINT 

PRINT 

PRINT 

PRINT 


... 

pi?'* 
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ORIGINAL PAGE IS 
OF POOR OTTaLTTY 


PRINT */'Prof il# (»lpna)«'/ALPHA 
PRINT */'T infinity **/Tinf 

PRINT */ 'Chas sal ra si s t . = */ RE $1 ST ANCE 
PRINT */'Nu*. points =*/NCPTS 

PRINT */'0ai Sate = * / OElSET A 

PRINT # 

PRINT */' Harkal' 

PRINT »/'csnstants NS = '/NS 

PRINT */' ME = */ME 

print Alp s'/HGHT 

PRINT », 'integration b ght * * /HGHT* 0EL5 ET A 

PRINT 

RETURN 

EN3 


SuERCLTINS SLILGFaTRIXCDELESTA/ALP) 




C 

C 

C 

C 

C 

C 


c 


c 

c 


* Tnis subroutine is used to determine the function 

* "G" in Van Iscrnd's notes. It «ill b e called by 

* the rain program and will r e t urn the matrix eh ich 

* mil c« inverted by tha Hankel program. 


IMPLICIT CCU5LE P R =CI SION ( A-«/0-Z) 
CCMMCN / MAIN/ C-oAR 


COMMON /CONST/ CI/PI 

CCMMCN /CONST 1 / irfcC/OMEGA 

CCMMCN /CONST!/ alPHA/OTCT 

COMMON /CON ST 3/ Z/ZREF/S 

CCMMCN /REGION/ R $1 / R S2/R II /R 12/ R 01 / RO 2 

COMMON /STATE/ S OE/ R E S 1ST A NCE 

CCMMCN /INTEG/ NS/M£,N9PTS 

C0M?LcX*1o GSAR(2276S)/3eTA 

C CMPL E X*1 6 GZ_Z/CZ_S/GZ_C/OZZ_0/GZZ/GZ 

C CMPlE X *1 6 ETA.Z/STA.S/ETA.O 


CGMPLEX*1 6 Ml _NZ/h2.NZ/Hl_N$/H2.hS 
CCMPLEX*16 Hi. NO/ M2.NC/ Mil. NO/ H21.N0 
C0MPLEX*1 6 CONST /R/E/CI/T1/T2/T3 
CQMPLEX*1 0 TAU/PS'l/ZIMP 
C0M?LEX»16 SRGZZ/SRGZS/SRGZ/SRGS 
INTEGER IREGION 

E = COEXP (* PI * Cl / 3.00) 

FR * OMEGA/ C2.D0*PI) 

RATIO = FR/RE5ISTANCE 

CO = 1 .00*9. C80C*RATI0**C-. 7500) 

X = -11.90C*RATI0**(-.73D0) 

ZIMP * DCMPLX < 00/X ) 

Q * T .00 

a S 0 . C 0002 *4 . CO *PI* 4 .47D0* (1 0 . O0**< SPLRE F/20.00) ) 

RLMOA » OMEGA/(SPEED*ALPHA) 

RLM0A23 = RLMCA ** (2.00/3.00) 
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************** ************** 

*• PROGRAM 8EGINS HERE ** 

CO 1C/I=1/NCPTS 

o£TA = OCKPIX CDFLCAT(I-1)/(ALP))*(DEl.BETA) 

CAUL G Z ALLl CZREF/3ETA/GZ.C/GZZ.0/ETA 0) 

CALL GZALL1 ( Z/ 3 ET A/ GZ. Z /G Z Z /ET A Z) 

CALL GZALL1 CS/EETA/GZ.S/GZZ/ETA.S) 

CALL HALL (ETA.C/H2.NC/H21 _NC/H1~N0/H11 _NC/ IREGO) 

CALL riALL2CETA.$/H2_K3/Hl_NS/IRSG$) 

CALL HALL2C5TA_Z/H2_NZ/ri1 .NZ/IREGZ) 

CALL REGION C £ E T A , I R c G I ON ) 

TAL = 1.CC - CSPEEC/CH5GA * CI/2 * ZIMP * GZZ.O/GZ.O) 

3 si = :ihp*ci*speec/chega 

-SI = P 3 1 * CC3.CC/2.C0 * RLM0A)**(2. 00/3.00)) * GZ.C 

SAGZ = COS CRT C G Z _ Z ) 

SRG3 = COSCRT < G Z _ S J 

CONST = C * EETA / Cl 2* C I * ( R L MC A2 3) *S RG Z *SRG S ) 

If CCAcSCri2.N0)) .EC. C.DO .ANC. CA8SCH21.N0) .EQ. 0.00)) THEN 
R = 0CMPLX(1 .035/0. CO 
SOTO 22 

ENG IF 

R * TAU « H2.NC + PSI * H21.N0 

R = -CTAO * "1.NC + 3 5 1 * H11.N0) / R 

SC TCC1/2/3/A/3/t/7/8)/IRESICN 


• RE SION 1 • 

GEAR (I) = CONST * CHI. NS ♦ h2_NS*R) * H2.NZ 

SOTO 10 

• REGION 2 • 

GoASCI) = CONST • CH1.NZ + H2_NZ*R) * H2.NS 

GOTO 10 

• REGION 2 * 

G3ASCI) s ~ (CONST /(E*E>) * CHl.NS/R * H2.NS) * H2.N2 
GOTO 10 

• REGION 4 • 

GsARCI) =-CCNST ,* CHl.NS/R ♦ H2.NS) • H1.NZ 
GOTO 10 

• REGION 5 * 

G3ARCI) =-CCNST • CH1.NZ/R ♦ H2.NZ) * Hl.NS 

GOTO 10 


** REGION 6 ** 

GSARCI) a CONST • (Hl.NS ♦ <<R-E>/R)*E • H2.NS) • H2.NI 
GOTO 10 

* REGION 7 * ~" 

G3ARCI) * CONST * CH1.NZ ♦ C(R-E)/R)*E * H2.N2) * H2.NS 

GOTO 10 



* 


REGION 8 * 



ortoooonon n o o r> f) o o n r» n n r> n 


ORIGINAL PAGE IS 
OF POOR QUALITY! 

G8 AR ( I ) « CONST * E * CH1.NZ/R ♦ H2_NZ) * H2.NS 


10 


CONTINUE 

RcTURN 

£N0 


SUBROUTINE G32ALL ( 2 / B E T A,G3 2C ) 

A************************************************************* 
* * 

* G3 2 ALL CALCULATES g3/2 FUNCTION * 

* * 

* * 

* mcaified so the oranch point of the LCG function * 

* is below the positive real axis. F*b 20* 86 * 

* modified again 1C/12/86 for mocified progn* * 


***x«*ttKw************#ft*****fc#*«t**+*******4*#****##******#** 

IMPLICIT DOUBLE PRECISION (A-H/C-Z) 

C0M?LEX*16 SET A, PHI/ S 0 R T 1 / $ CR T 2 / F 00/ G 3 2C / MOOLOG 

C0MPLcX*16 S1/A/2/S/CI 

COMMON /CONST/ Cl/rl 

COMMON /CCNST1/ SPEEC/OMEGA 

COMMON /CONST!/ ALPHA/OTCT 

A A = 1.00 * ALPHA*! + DTCT 

3 = 1.00 - ( (S?EE0*BSTA/CMEGA) * ( S PE EC * B E T A /OMEG A) ) 

SI = S i 3 T 1 (EETA/I) 

S = SOSTKcETA) 

Phi = SI /(S*CSQPT < A A ) ) 

IFCGTCT .EC. ,C .OR. CD AB $ < 1 . - t>Hl ) .LT. 1.0-8) THEN 

“03“ 0 . vj 

E kSE 

F00 = .5 * CCLCGC (1.+PHI)/(1.-PHI> ) 

ENDIF 

G32C = - OSCRT(AA) * SI + CTOT • POO/S 
RETURN 
c NO 


BEGIN OF FUNCTIONS USEC ABCVE. 


SCRT1 AND SCRT2 ARE FUNCTIONS TC CALCULATE SORT (BETA) GIVEN * 

CESISED BRANCH CUTS AND 3IRECTI0N ♦IMAGINARY CR -IMAGINARY. * 

ft*/*/************************************************************* 


FUNCTION SQRT1 (SETA/Z) 

IMPLICIT OOUBLE PRECISION CA-H/O-Z) 

CCMPLEX*1 6 3ETA/AA/SCRT1/CI 
COMMON /CONST/ CI/PI 
COMMON /CCNST1/SPEED/CMEGA 
COMMON /CONST2/ALFMA/CTOT 

AA * 1.-((CSPEEC/OMEGA)*BETA)*<(SPEEO/OMEGA)*BETA)) 

83 * 1. ♦ ALPHA*Z ♦ OTOT 

SCRT1 *CDSCRT (AA*BB-DT0T) 

RETURN 

END' " 
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FUNCTION SQRT2 (BETA) 

IMPLICIT OOU3LE PRECISION (A-H,0-Z> 

CCMPLEX*1 6 SETA/A*, SORTS, Cl 
COMMON /CONST/ Cl, PI 

COMMON /C0NST1 / SPEEC, OMEGA 

COMMON /CONSTS/ ALPHA, CTCT 

AA = 1.-(USPEcC/CMEGA)*9ETA)*((SPEE0/0MEGA)*86TA)> 
SORTS = COSORT (AA) 
kETURN 
cNO 


SuSROUTINE GALL( Z,5STA,G) 

*****,**»»,. ******•********, ******** ************** 

GALL EVALUATES ThE g FUNCTION * 

IMPLICIT OCUSLE PRECISION (A-h,0“Z) 

COMPLEX *16 C-, SETA, 532,32, Cl, gT 
COMMON /CONST/ Cl, PI 
COMMON /CCNST1/ SPEEC, OMEGA 
COMMON /C0NST2 / ALPHA, OTCT 
CALL GZ2AILCZ, SETA, 532) 

PHI = OATAN2 (OIMAG(G3S),OREAL(G32) ) 
if (phi .gt. 3 . d C ) phi = oM-2*pi 

G = (Ccdabs(932))**(2.dO/3.dO)) * cd**p < <phi*2.dO/3.dC)*ci) 
si =3 

IF (Phi .Lc. -PI/Z.DCT THEN 

G = G * CCEXP(2.CC/3.30*PI*CI> 

ELSE 

G = 3 * CDEXPC4.DC/3.00*PI*CI) 

ENC IF 
RETURN 

END ' ' ~ " . — — - 


C 

C 

c 

c 


SUBROUTINE GZALLTCZ, BETA, GZ,GZZ,EW) ~ 

********************** ******************** ************** 

GZALL1 CALCULATES ALL THE PARTIAL DERIVATIVES * 

CF THE g FUNCTION. *' 

********************************************************** 
IMPLICIT OOUeLE F REC I S I ON ( A-H,0-Z ) 

COMPLEX *16 3cTA,GZ,GZZ,G,GB,GBS,SQRTT~, EN,CX,OUM 
COMMON /CONST/ Cl, PI 
COMMON /CONST1/SFEEO, OMEGA 

COMMON /CONST2/ALPHA,CTOT 

CALL GALL(Z,BETA,G) 

SI « 1 • DO 

A * 1.*alPHA*Z*OTOT - 

GZ « 2.00 * SI* ALPHA * SQRT1(BETA,Z) / <3 . *C DSQ R T ( 6* A ) ) 
C' • 2.*ALFHA**3.0C«DT0T/<9*A**2.00» 

GZZ - C/ (GZ»6J-*.3*GZ**2.D0/C * 

T * <3.*0NEGA/(2.*ALPHA»SPEE0))**<2.OCA3.0O> 



i N = T*G 

RETURN 

END 


ORIGINAL PAGE IS 
OF POOR QUALITY 


SU 3RCU T IN E HALL (Z/H2/M21 /HI / h 1 1 /1REG) 

r * * * * 1 

nA u L U$=3 SUBROUTINE CG3AIR TC CALCULATE 1/3 CROER 
h ANNEc FUNCTIONS FROM AIRY FUNCTIONS. 


C»#*.»»**»»»*»*»*»»»»**»».*****»< 


IMPLICIT C0U3LE PRECISION (A-H/O-Z) 

COMPLEX* 1c I/AZ/eI/Al?/3IF,X,K$/H1,H2/Hl1/H21 

uCMPlcX*1c ARC-/CI 

COMMON /CONST/ CI/PI 

ARC = CCN?LX(C.CG/-PI/6.00) 

,< = C12.CG>**C1 .00/6 .OC)«COEXP(ARG) 


xS = OCCNJGCX) 

CALL CCe A IP (-I/AI/el/AIP/’SIP/IREG) 

nl = K»(AI-CI*cI) 

nZ - xS*(4I+CI*cI) 

nil * -X»(AI?-CI*EI°) 

nil = -XS*(AIP+CI*sIP) 


RETURN 


END 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUo SOU TINE CG6AIR<Z/AI/3I/AIP/8IP/IReC-> 
.*»*»»**•«***«*********< 


CALCULATE AIRY FUNCTIONS FC R CCMPLSX*16 ARGUMENT * 

REF. .-AN03C0X OF MATHEMATICAL FLNCTICKS/ A3RAMOWITZ ANO STEGUN. * 
ENTRY: * 

CALCULATE ARGUMENT <Z) ANC ASSCLUTE VALUE(Z) * 

IF 111 LT 6 * 

Then USE cCS. 1C. 4. 2 ThRli 10. 4. 5 FOR Ar/31/AIP/eiP * 

10 Et-cE IF ARGCZ) LT PI/3 * 

THEN CALCULATE ZETACZ) * 

USE SIS. 10.4.59/ 10.4.61/ 10.4.63/ 10.4.66 FOR AI/ 8 X/AIP/ 6 IP* 

20 ELSE CALCULATE ZETA(-Z) * 

USE ECS. 10.4.60/ 1C. a. 62/ 10.4.64/ 10.4.67 FOR A I / 81/ AIP/ 8 X P* 
ENQIF * 

ENOIP * 

EXIT ' * 

END * 


IMPLICIT C0U8LE PRECISION (A-H/O-Z) 

COMMON /CONST/ CI/PI 

C0MPl£X* 1 6 Z/AI/3I/AIP/BIP/ZETA/CZETA/Z14/SUM1/SUM2/SUM3/SUM4/ 
1 ZETA?/FACT1/FACT2/SN/CS/FTERM/FPTERM/GTERM/GFTERM/P/FP/G/GP/Z3 
COMPL EX*1 6 VZETA/VZETAP/CI 
OX McNSION C (21 ) /O <21 ) 

CATA C1/C2/PIRT/P 14/. 355028053900/. 25881 9403800/ 1.772453851 00/ 
♦ .735398163500/ 

? IRT s GSQRT (P I ) 

OATA C/1 .00/.C694444444444440C/ 

+ .C 37 1 3343765 43 21 OC/. 03 799305 91 2 7 SCO 00/ 

1 .C5764919041Z6693C/.11 6099C64C255100# 
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♦ .291591 3992307400/ .87766 69 6 9 5099 8C0/ 

2 3. 07 94 5 30301 73 100/ 12. 341 573 3 32 3 4 5C0> 

+ 55. 62278536591 400/ 278. 465C8077759C0/ 

3 1533. 16943201 27 00/92 07. 2 06 5997 258C0/ 

♦ 55892. 51356587500/419524.87511653C0/ 

4 31 48257. 41 78666C0/ 251 58915. 871601 CO/ 

♦ 21 4238036. 9636600 /I 929 375 5 49,1 8 23CO/ 

5 18335766937.88900/ ' 

3ATA C/1.00/ 

♦ ”. 097222222222221 CQ/-.C4388503C86419700/”. 042 462 83 078989400/ 

1 -. 06266 21 6 34 92 031C0/-. 1241 0585 60 272700/-. 3082 53764901 0700/ 

2 -. 9204799 y 241 29 1 CC/- 3. 2 1049 358 46 48 500/-1 2. 80729308073500/ 

3 -57. 5C3 3035 1 39 1 1 0C/-287. 033 2371 092000/-1 576.357303337000/ 

4 -944s. 3 5 <.£2 309 530C/-&1 3 3 5. 70666 334700/ -4 289 52.4004000400/ 

5 -3214536. 521 4C060C/-256979C8.3839C900/-218293420. 8321400/ 
c -1 963 523788. 99C90C/-1 8643931088.10500/ 

A6S2 * A s S ( 2 ) 

IF(ASSZ.EC.Q) GO TO 3 

IF(AaS(0IMAG(2)).L6.1.C-12. AND. OREALCD.LT. 0.00) GO TO 5 
ARGZ » CATAN2(DIM6G(Z)/0REAI(Z>) 

GO TO 4 

3 ARGZ > 0 . DO 

GC TO 4 - 

5 ARGZ * o i 

4 CONTINUE 

IrCASSZ.GT .6. CO) CO TC 10 "" 

ir *g*1 

C ASCsNCING SERIES 

C 80S. 10.4.2/1C.4.3 

15 CONTINUE 

13 a i**3 

FTERM a t.DO ' 

FPT5RM* Z»I/2.00 

GTERM = Z 

GPTERM* 1.00“ - - 

GLIM * 1.D-13*ABSl 

F a FTERM 

rP ■' * FPTSRM - - — 

G » GTERM 

GP a GPTERM 

KKKT a ICO -f 'ADJUST RKICT TO INSURE CONVERGENCE IF NECESSART 
00 1 1*1 /KKKT 
13 » 3*1 

FTERM- « FTSRM*r3/T ( I3-r.DC)*I31 * ' ~ 

FPTERM=FPTERM*Z3/ (13* (1 3*2. 00) ) 

GTERM * GTSRNAZS/flSaCISM.OO)) 

GPTERFaGPTERM*Z3/C(IJ-2. 00*13) 

F * F+FTERM 

FP * FP+FPTERM 

G * G+CTERM ' “ 

GP * GP+GPTERM 

IFCC0A8S (GTERM) .LE. GLIM) GO TO 2 

1 CCNTINU E — — 

PRINT 6000/ 1 

6000 FORMAT (/' Z**2fU.5/* ERROR IN CGBAtR/ NONCONVERGENCE*) 

2 AI — m— C1*F-C2*C 

AIP » C1*fP-C2*GP 

81 * 1 . 732C508C800*(C1*P*C2*6) 

8IF * 1.732C508C800*(CT*rFK2*CF> ~ 

GO TO 9999 - ,j 4 

■■ 4 % -- ■ 
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ASYMPTOTIC EXPANSIONS FOR ✓ Z / LARGE 
SIGN * 1.00 
SUMl * 0.00 
SUM 2 * 0.00 


SU M3 = 0.00 
SUM4 = 0.00 

IF ( A 3 $ ( A R G Z ) .GE. PI/3. DO) GO TO 20 
l r a 3 = 2 

/ ARG ( Z ) / Lc PI/3 

cCS. 1C. 4 . 59, 10.4.61, 10.4.63, 1C. 4. 66 



2c T A 

CC 11 1= 
K - 

2 c T A P = 
SC Ml 
SUM 2 

S L' M 3 - 

Surt4 = 
- 

Con tirua 
2 U 

r A CT 1 * 

FACT 2 » 

A I 

Ml ? 

r ACT 1 = 
FACT2 = 
o I - 

6 1 P 


CZcTA CA6S2/ARG2) 

1,12 

1-1 * 

:eta**k 

$ JM1+SIGN*C(I)/Z5TAP 
SUM2 4 $ISN*C(I>/ZETAP 
SUM3 4 C(I)/2ETAP 
SuM4 + 0C)/2cTAP 
-SIGN 

A5S2 **.2500 * DCPPLX(COS<ARGZ/4.00>*SIN<ARGZ/4.CQ>> 
.500* COcXPC-ZcTA)/CPIRT*ZU> 

.500* C0EXP(-ZETA)*Z14/PIRT 
FA CT1*$UP1 
-FACT2*SLM2 
C3SX?(ZcTA)/(PIRT*Z14> 

COEXP (ZcTA)*Z14/FIPT 
FACT1 * SUP 3 
F AC T 2 * SUP 4 


SC TO 9999 


/ARGCZ)/ GT PI/3 NOTE CHANGE A 8 CVE 
cwS* 10.4.60, 10.4.62, 1C.4.64, 1C.4.67 


CONTINUE 
ir eg = 3 

A RGZ = OATANZC-OIMAGCZJ,-DREALCm 

ZcTA = CZcT A<A 6 SZ,ARGZ) 


VZETA = 1 .OC/ZETA 

LLL = 10 

00 21 1=1, LLL 
K 2 = ( 1 - 1 ) *2 

J = K2+1 

VZcTAP = VZETA**tL2 

SUMl = SUM1+SIGN*CCJ)*VZETAP 

SUM 2 = SUH2+(SIGN*C(J+1)*VZETAP*VZETAJ 

S UM3 = SUM3 + (SIGN*C(J)*VZETAP) 

SUM4 = SUM4+(SIGN*C(J+1)*VZETAP*VZETA) 

SIGN = -SIGN 

Continue 


ZU 

FACT1 

FACT2 

SN 

CS 

AI 

AIP 

ai 

6 IP 


= ABSZ**.250Q*DCMPLX(C0S(IRGZ/4.0C)/SIN<ARGZ/4.D0)) 
= 1 .OC/ (PIRT*Z14) 

= Z14/PIRT 
= SINCZETA+PI/4.Q0) 

= COS(ZETA+PI/4.DC) 

= FACT1 *<SN*SUM1-CS*$l!M2) 

= -FACT2*(CS*SUM3+SN*SUM4) 

* FACT1*<C5*$UM14SN*$UH2) 

* FACT2* (SN*SUM3"CS*SUM4> 
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C 

9999 continue 
RETURN 
ENO 


ORIGINAL PAGE IS 
OF POOR QUALITY 


BEGIN OF FUNCTIONS USEC ABOVE 


FUNCTION CZETA <ABSZ*ARGZ> 

IMPLICIT COUalE PRECISION CA-H,0-U 
CCMPLEX*14 CZETA 
ARC = ARSZ*1.SOC 

CZETA = (ABSZ**1.5CG)*CCMPIX(C0S(ARG) / SIN ( A RG) ) * . 66666666E 6 666 700 

RETURN 

ENO 


ENC OF FUNCTIONS LSEO ABOVE. 


SUBROUTINE h AL L 2 ( Z/H2, hi , IR EG) 

C »••*»»««»*****»******»*********»**************#******#****** 

C nini USES SUBROUTINE CG5AIR TC CALCULATE 1/3 OROER * 

c t-ANNSL functions from airy functions, nct the * 

C DERIVATIVES as hall cces. * 

c»*»»*»»«»* »**»«*»»»»«*«***»«****»*»*»*##**»*****#******#**##* 

IMPLICIT C0U5LE PRECISION (A-M/C-Z) 

COMMON /CONST/ CIrPI 

CCMPlEX* 1 6 Z, AI,aI,AIP,3IP/K/l(S*H1/H2,Hl1,H21 
CCMPLcX*1 6 ARGsCI 
CCMPLEX*16 beta 
ARG * 0CKPLX(0.0C*-PI/6.CC> 

X * <12.0Q)**(1.00/6.DC>*C0EXP<ARG) 

NS « CCONJG(K) 

CALL CGSAIR(~ZsAIs3I/AIP/BIP/IREG) 

nl * K*(AI-CI*3I) 

m2 * KS*<AI+CI*EIJ 

RETURN 

ENO 

C " 

C 

C 

C • -■ ' ■ 

SUBROUTINE HANKEL 

C *•*•»*• *.*****•••***••« •*»***»**»»•**** *»** a************************ 
C SUBROUTINE TO ACCURATELY CC THE HANKEL TRANSFORM CF THE SOUND 
C PRESSURE LEVEL. 

C F ( NP) S GBAR (NP ) MUST EE SAMPLEO AT NF POINTS WITH K«(N-1,ALP) 

C ALP REPRESENTS THE OISTANCE ABOVE THE REAL AXIS THE FUNCTION MILL 

C EE INTEGRATED. 

C NS IS A PARAMETER REPRESENTING AOOZTION OF AN ANALYTICAL FUNCTION 
C TO FCNPJ r"' — - ' 

C M IS THE NUMBER CF TERMS USEO TO APPROXIMATE F (NF) TO INFINITY 
C*»»* * ******* *********** *****••»•••**»••« •••*••*«****•* *************** 

IMPLICIT DOUBLE PRECISION '7 — ~ 

COMMON /CONST/ CPPI/PZ 
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C0MMCN/AF82IN/ALF/0EL6STA 

C0M?L£X*16 FU274 8) / CF/CARG/SUM/FNP/CMPI/OI 

ALP * -ALP ! naca**ary to aatcN attanberough'a definition of alp 

HP * HI 


C 

c 


OELK = 0EL3ETA 

SLe TRACT THE ANALYTICAL FUNCTION IF NS > ZERO 
ACJoSTlNG THE SUBTRACTION MULTIPLIER CF 
IF(NS.LE.O) GOTO 11 
CF * OCMPLX (0.CG/0.C0) 

IF (ALP.Ei.C.C) then 
CF = GFLOAT (N» > / C PIC AT (NS > 

Cr * C F *r ( 2 ) 

EnD if 

IF ( A L P . N E . 0 . C ) TrEN 

CF = CMP:*CFLQAT(NF)*F(1)/(CFL0AT<NS)*ALP> 
END IF 

iUsTRACT the ANALYTICAL flncticn if ks>o 
CC 1C,I=1,NP 

Cl - OCMPLX (CFLCATC-1 ), (-ALP)) 

CARG = CFLCaT(N$)=(-C1)/CFL0AT(NF) 

S (I> * F(i)-CF*(1.C:-C0EXP(CARG)> 

1 J CONTINUE 

11 CONTINUE 

Ir(ALP.SC.O.C) F(1)=CCMPLX(0.00/C.DO) 

FnP = r ( N° ) 


12 

C AGC 


15 

2G 

C DC 


C ADC 


25 


CC 1 2,C*2,N? 

01 = DC MR LX (O-LCAT (1-1 ) , (-ALP)) 

.-<!)= FC)/(C0SIRT(C1)) 

CONTINUE 

IFCALP.N5.0.C) =(1)=F(1)/CCSCRT(C-CMPI)*(ALP)) 
TeR.yj TC INFINITY IF HE>0 
IF(ME.LT.I) GOTO 20 
CC 15/1=1, NP 

01 = CCHPLX(DFLCAT(I-1>/(-ALP>) 

CF = 01 /C c LOAT (NP) 

CALL IcTA (NP/Mc/CF/SUM) 

F ( I ) = FCI)*FNP»5LM 
CONTINUE 
CONTINUE 
The FFT 

CALL FORK (NP/F/ 1 ) 

ALTERNATE TERMS TO GIVE NP/2 SAMPLES TRANSFORMED 
Cf = OSLK*CFLOAT(NP)*(CCSCRT(-CMPI))f (2.00*FI> 
CO 25/I=2/NP/2 

A = DEXP(CFLOAT(I-1)*(ALP)*2.CO*PI/DFLOAT(NP)) 
r ( I) = A*F(I)+CMFI*F(NP-I+2)/A 
F(I)= F(I)*CF/0S0RT(CFLCAT(I-1)) 

CONTINUE 

RETURN 

END 


SUBROUTINE FCRKCLX/CX/SIGNI) 

•••••••I************************************************************ 

A FAST FFT GIVEN SY J.F. CLAER80UT/ "FUNDAMENTALS OF GEOPHYSICAL * 
CATA PROCESSING" * 
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IMPLICIT DOUBLE PRECISION CA-H,0-Z> 

COMMON /CONST/ CPPI/PI 

COMPLEX* 1 6 CX(LX), CARS, CW, CTEMP, CI,0UM1,CMPI 

INTEGER SIGNI ' 

J * 1 

SC 3 0$QRT<1.00/CFL0AT<LX>) 

00 3C,I*1/LX 

IF(I.GT.J) GOTO 10 
CTEMP * CX(J)*SC 

cxtjj * cxci)*sc ' • - • 

CX<I) » CTEMP 
10 M»LX/2 

20 IF(J.LE.M) GOTO 3C 

J * J-M 

M = M/2 

IFCM.GE.1) GOTO 20 
30 J = J* M 

L - 1 

40, 1 3T5 P = 2*L 

CO 5G,M=1,L 

CASG = CMPI*PI»0bLi:(SIGNI)*C3LE((M-in/0BLE(L> 
C*=CCEXPCCARG) 

00 5C,I=M,LX,I$T6P 
CTEMP* CN*CX<I+L> 

CXC:+L)=CX(I)-CTEMP 
50 CX CI)=CX C I) ♦CTEMP 

L = I 5 T e P 

if(l.lt.lx) gcto ;c 

R ETuSN 
END 


SU5RCUTINE ZETA ( X P, M/ A/ SUM ) 

•**«******************«********************************* 
SUBROUTINE TC AQC THE NECESSARY TERMS TO 
EXPRESS INVERTIBLE FUNCTION TO INFINITY. 
wlLL USE DOUBLE PRECISION. 

SUM=SUM OF 1 /(NP*.5)*1/CCJ+A)».5> FOR J*T TC 
INFINITY MINUS SOME CONSTANT WFICH IS 
INOEPENCcNT OF A. 

a******************************** *********************** 
IMPLICIT 00U6LE PRECISION (A-H,0-Z> 

C0MPLEX*16 A,SUM/C2 

02 * OCMPLX (DFLOATOO ,0. CO) 

SUM * 1 .00/ <M*A) 

SUM * 2.00*(0SQRTC0FL0AT<M))-1.C0/C0SQRT($UN)) 

S -0.3*C0S0RT(SUM)*(1.C+SL'H*(1 .0/12.0-SUM* SUM/1 92. OJ) 

CO 1 C, J*1 ,M 

SUM * SUM+1.00/CGSQRT(J+A) 

10 CONTINUE" "" - 

SUM * SUM/OSQRT(CFLCAT(NP)) 

RETURN 


ononoonnr>oo o o o r> ooonoooo 


SueRCUTlNE PRINT4LL(C6L3'Tfl/N) 



* 


SUBROUTINE TO PRINT THE RESULTS OF TFE 
HANKEL TRANSFORM TO THE FILE FORC44.0AT 



IMPLICIT OOUSLE PRECISION <A-H,0-Z> 

COMMON /CONST/ CI/PI 
CCMMCN /MAIN/ F 
CCMPL£X*16 C <32763)/CI 
CO 1CC/I=1/N/2 

RAC = 2.3G*P> CFLOAT(I-I) / ( C FL 04 T <N ) *0 ELS ET A ) 

IF C R AO .ST. Q.O) THEN 

OcCIBLE = 2C.CC * CL0G1C(CCA3S(F(I))) 

RA ZZ = uLOGICCRAO) 

wRITE ( •* 5 ✓ 9 CC ) RAD/RA02/:ECI3LE 
E N 0 Ip 

too continue 

*uu FORMAT (5X/5G15.7) 

RETURN 

end 


SUBROUTINE REGICN_ C INC C A N G/ R T 1 , R T2 ) 


(*+****++< 


Inis sucrcutm# determines inhere region changes 
in the gI/2 function take piece. The routine eill 
be csllec by the aain prcgrae and eill return 
tne varieties RSI / R S2/ RZ1 / R Z2 / RQ1 and R02 which 
are the values of beta where region changes occur 


*«**•**»* 


5 


1C 


implicit acueus precision ca-h,o-z> 

CCMMCN /CONST/ CI/PI 

CCMMCN /CONST 1 / SPEEC/CMEGA 

CCMMCN / CONST 2/ 4LPHA/DTCT 

CCMMCN /CQNST3/ 2/IREF/S 

COMMON /REGION/ RS1/RS2/RZ1/RZ2/R01/R02 

COMMON / AF32IN/ FGHT/CELEETA 

COMMON /INTEG/ NS/ME/NOPTS 

C0MPLEX*1 6 3/C I 

DIMENSION RS(3)/RZ(3)/RQC3) 

print */'first rcot is'/RTI 

print «/'s»ccnd root '/RT2 

HLIM * RT2 * 1.C100 

SLIM * REAL(INT(.95C0*RT1/DEIBETA>) * OELBETA 
M * HGHT *0EL8£T A 

i * o " 

J = C 

K *0 

B = OCMPLX(BLIM/H) 

CALL G32REGI0N(S/8/I/ANG/PS1) 
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55 


t 


1 


CALC G32R£GICN(Z,3r J,ANG/PZ1 > 
CALL G32REGICN(0/3sKsANG«P01) 

R S ( I ) » REAL(d) 

R Z ( J ) * REAL (3 ) 

RGOO * SEAL(a) ' 

3 = 3 ♦ CEL3ETA 

IP UDRE6LC3) .Lc. ML I*' ) .ANC. 
print 

IP (ItJ + K .LT. 9 ) P R I N T *,* 
R$1=RS<1) 

RS2*R$(2) 

RZ1*RZ(1 ) 


original- pags i3 

OF POOR QUALTra 


.LT, 9>) GOTO 10 

REGION BOUNDARY NOT 0£TECTE0# # 


RZ2=RZ(2) 

aoi=RC(l ) 

R02*RC(2) 

NT * :i.nT((NIN(RS 2^R22/RC2)-R01)/D6LBETA) 
if (nt • £t . C) than 

PRINT */'#st, Miri/nu* nuTb«r of points in rsgion 2 is*#NT 
ana if 

print * * * R 31 / R SI '✓RS1#RS2 
print */'RZ1,RZ2 ',RZ1,RZ2 

print * * ' RG 1 / RC2 # /RG1/RC2 
R c T 0 R N 
c NO 


i 

1 


SUBROUTINE G22RcGION( Z / 3 ✓ L / ANG / Pi ) 

*******x «*******-*»**************4 

* Subroutine tc C2T=rkinc g 32 function and where region changes occur * 


[***★***< 


IMPLICIT CCU3LE = S E Cl S ION ( A 0-Z) 
COMMON /CONST/ CI/PI 
C0MPLZX*16 3,G32,CI 
CALL G32ALl(Z,S/G32) 

P 2 = 0ATAN2CQIMAGCG32)»CREALCC-32)) 

IF (P2 .GT. C.OO) P2*P 2-2*PI 


IF (L .EC. 0) 

THEN 



L = 1 

GOTO 1 




END IF 




IF <<P1 .LT. ANG 

.ANC. 

P 2 ,GE. ANG) .OR. 


(PI .Gc. ANG 

« A NO * 

P2 .LT. ANG)) 

L*L*1 


PI * ?2 
RETURN 


SU3RCUTINE REGION C BETA / IR EGICN) 


* This suoroutin* Bill dtttrsins ihieh region is * 

» currently bsing addrtjssd By tBs progra. * 

******•»•*******«****»•••»*•***»*•»»***»**•»*•••***•*•*»**••***** 
IMPLICIT QOUBLE F RECIS ION < A-H, 0-Z) 

COMMON /CONST3/ Z/ZREF/S 

COMMON /REGION/ R SI ,RS2r R Z 1 / RZ 2*R01 ,R02 
COMMON / APB2IN/ PGHT/CELBETA 

C0MPLEX*16 BET* ~ 

3 • OR EAL(BETA)-OEL6 EtA/5.00 
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IF (Z .GT. S) then 









IREGION 

* 

t 

IF 

((a 

.GT. 

ROD 

.ANC. 

(e .is. 

RC2)> 

IREGION 

* 

6 

IF 

((6 

.GT. 

RSD 

• amc • 

CB • LE « 

R S 2 ) ) 

IREGION 

* 

3 

IF 

((a 

.GT. 

R 2D 

.AMO. 

(8 .LE. 

RZ2)> 

IREGION 

* 

4 


else 






IREGION 

X 

2 

IF 

cca 

.GT. 

ROD 

.ANC. 

(8 .LE. 

RC2) ) 

IREGION 

X 

7 

IF 

( ( B 

.GT. 

R 2D 

• ANC* 

(3 .LE. 

R 2 2 ) ) 

IREGION 

X 

8 

IF 

((8 

.GT. 

RSI) 

• ANC. 

(3 .LE. 

R S 2 ) ) 

IREGION 

X 

5 


EMC IF 
RETURN 
END 
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